Environmental and local habitat variables as predictors of trophic interactions in subtidal rocky reefs along the SE Pacific coast

Clear the environment
Load packages
Part 1: Importing and examining data
- 1a) Import data
- 1b) Inspect structure of data frames
Part 2: Manipulating datasets
Part 3: Latitudinal patterns in explanatory variables
- 3a): Explore latitudinal relationships
Part 4: Latitudinal patterns in consumption
- 4a): Predation vs latitude
- 4b): Herbivory vs latitude
Part 5: Predation intensity: Effect of site-specific variables
Part 6: Herbivory intensity: Effect of site specific variables
- 6a) Model selection
Part 7: Visualisation of best fitted predictors
Part 8: Validate sea water parameters
Part 9: Principal component analysis (PCA) of predictor variables
Part 10: Size estimations summary table
Part 11: R Session Information

Clear the environment

rm(list = ls())

knitr::opts_chunk$set(warnings=FALSE)

Load packages

Before loading the packages install them using install()

library(tidyverse)
library(summarytools)
library(dplyr)
library(readr)
library(vegan)
library(patchwork)
library(pander)
library(ggeffects)
library(performance)
library(knitr)
library(car)
library(corrplot)
library(MASS)
library(MuMIn)
library(geosphere)
library(mgcv)
library(glmmTMB)
library(lme4)
library(factoextra)
library(fitdistrplus)

#If Issues installing packages
options("install.lock"=FALSE)

Part 1: Importing and examining data

In this part the data frames were inspected in order to check the variable names and types. In some situations the data set was slightly modified for further manipulation of them.

1a) Import data

Environmental1 <- read.csv("Data/Environmental_variables.csv")
Metadata <- read.csv("Data/Metadata_lat2018-2019.csv")
Herbivory <- read.csv("Data/Herbivory_latitudinal2018-2019.csv")
Benthic <- read.csv("Data/Percent_cover.csv")
Phlorotannins <- read.csv("Data/Phlorotannins.csv")
Predation <- read.csv("Data/Predation_latitudinal2018-2019.csv")
Size <- read.csv("Data/Size_estimations.csv")
Transects <- read.csv("Data/Transect2018-2019.csv")
Sal_data <-read.csv("Data/extracted_sal_data.csv")
temp_data <- read.csv("Data/extracted_temp_data.csv")
o2_data <- read.csv("Data/extracted_o2_data.csv")

1b) Inspect structure of data frames

Site metadata

view(dfSummary(Metadata))

#Rename temperature variables to avoid confusion with replicate measures in Environmental data
Metadata <- rename(Metadata, Temperature_site = Temperature_C)

Environmental data

view(dfSummary(Environmental1))

Environmental <- Environmental1 %>%
  rename(D_oxygen_percent = X., D_oxygen_mgL = D_oxygen) %>% #Rename oxygen percent variable
  dplyr::select(-Date,-Campaign,-Latitude) #Remove variables that are in Metadata

Benthic cover

view(dfSummary(Benthic))
Benthic <- Benthic %>%
  dplyr::select(-Date) #Remove variables that are in Metadata

Herbivory

view(dfSummary(Herbivory))
Herbivory <- Herbivory %>%
  dplyr::select(-Date,-Campaign,-Season) %>%# Remove variables that are in Metadata
  mutate(Consumed_proportion = (Initial_length_mm-Final_Length_mm)/Initial_length_mm)

It was necessary to create major categories for Bait

Herbivory <- mutate(Herbivory, Bait2 = case_when(Bait == "Lessonia_spicata" ~ "Lessonia sp.",
                                                 Bait == "Lessonia_berteroana" ~ "Lessonia sp."))

Predation

view(dfSummary(Predation))
Predation <- Predation %>%
  dplyr::select(-Date,-Campaign,-Season)  %>%# Remove variables that are in Metadata
mutate(Consumed_proportion = (Initial_Individuals-Remaining_individuals)/Initial_Individuals, Consumed_individuals = Initial_Individuals-Remaining_individuals) #Add Consumed_proportion column

Create major categories for Bait (Squidpops vs Crabs) and Assay time (2 hours and 24 hours)

Predation <- mutate(Predation, Bait2 = case_when(Bait == "Squidpops" ~ "Squidpops",
                                                 Bait != "Squidpops" ~ "Crabs"), Assay_time2 = case_when(Assay_time == "24" ~ "24",
                                                 Bait != "24" ~ "2"))

Phlorotannins

view(dfSummary(Phlorotannins))
Phlorotannins <- Phlorotannins %>%
  dplyr::select(-Date,-Latitude) # Remove variables that are in metadata

Size

view(dfSummary(Size))
Size <- Size %>%
  dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata

Transects

view(dfSummary(Transects))
Transects <- Transects %>%
  dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata

Part 2: Manipulating datasets

In this section data sets are manipulated for analyses. New data sets with means and combined with metadata and consumption are created.

2a) Environmental data

Join Environmental data with metadata

Environmental_metadata <- Environmental %>%
    left_join(Metadata, by = "Site")

Join Environmental data with consumption

First summarise environmental data per site

Environmental_2 <- Environmental_metadata %>%
  group_by(Site, Latitude) %>%
  dplyr::select(-D_oxygen_percent, -Conductivity_Ms, -C_f, -C_i, -H_Visibility) #Keep only variables that will be used in the analysis

Environmental_2$Day_length <- daylength(Environmental_2$Latitude, Environmental_2$Date)


Environmental_sum_long <- Environmental_2 %>%
  dplyr::select(Site, Latitude, Temperature, D_oxygen_mgL, Salinity_ppt, V_Visibility, Day_length) %>%
  gather(Env_Variable, Value, Temperature:Day_length) %>%
  group_by(Site, Latitude, Env_Variable)%>%
  summarise(Mean_value = mean(Value, na.rm=T),n = n(), sd = sd(Value,na.rm=T), se = sd/sqrt(n)) 

Environmental_sum <- Environmental_sum_long %>%
  dplyr::select(Site, Env_Variable, Mean_value) %>%
pivot_wider(names_from = Env_Variable, values_from = Mean_value,
              values_fill = 0)

Then, add environmental data to consumption
Predation vs environmental data sets

Predation_env <- Predation %>%
    left_join(Environmental_sum, by = "Site")

Herbivory vs environmental data sets

Herbivory_env <- Herbivory %>%
  left_join(Environmental_sum, by = "Site")

2b): Benthic data

Join benthic data with metadata

Benthic_metadata <- Benthic  %>%
    left_join(Metadata, by = "Site")

Benthic_metadata_wide <- Benthic_metadata %>%
  group_by(Site, Latitude, Type) %>%
  pivot_wider(names_from = Type, values_from = Coverage,
              values_fill = 0)

Calculate mean, sd, and se by group

Benthic_mean <- Benthic_metadata %>%
  group_by(Site, Latitude, Type)%>%
  summarise(Mean_cover = mean(Coverage),n = n(), sd = sd(Coverage), se = sd/sqrt(n))

Join benthic cover data with consumption

First, summarise benthic cover per site and change it to the wide format

Benthic_sum_wide <- Benthic %>%
  group_by(Site, Type) %>%
  summarise(cover_mean = mean(Coverage)) %>%
  pivot_wider(names_from = Type, values_from = cover_mean,
              values_fill = 0)
Benthic_sum_wide_var <- dplyr::select(Benthic_sum_wide, Site, Kelp, Bare_rock, Understorey_algae)

Then, add benthic data to consumption
Predation_env + Benthic summarised data

Predation_2 <- Predation_env %>%
    left_join(Benthic_sum_wide_var, by = "Site")

Herbivory_env + Benthic summarised data

Herbivory_2 <- Herbivory_env %>%
  left_join(Benthic_sum_wide_var, by = "Site")

2c): Fish length

Join Size data with metadata

Size_metadata <- Size %>%
left_join(Metadata, by = c("Site"))

Calculate mean, sd and se per group

Size_mean_sp<- Size_metadata %>%
  group_by(Site, Latitude, Species, Type)%>%
  summarise(mean_Size_sp = mean(Size_cm))
Size_mean <- Size_mean_sp %>%
   group_by(Site, Latitude, Type)%>%
  summarise(mean_Size = mean(mean_Size_sp),n = n(), sd = sd(mean_Size_sp), se = sd/sqrt(n))

Calculate mean of fish size and add it to consumption

Size_mean_wide <- Size_mean_sp %>%
  filter(Type == "Fish") %>%
  group_by(Site, Type)%>%
  summarise(mean_Size = mean(mean_Size_sp)) %>%
  pivot_wider(names_from = Type, values_from = mean_Size,
              values_fill = 0) %>%
  rename(Fish_size = Fish)

Add fish length data to consumption
Predation_env + fish size data

Predation_3 <- Predation_2 %>%
    left_join(Size_mean_wide, by = "Site")

2d) Species list (abundances and richness)

Add metadata to Transects data set

#Join transects with metadata
Transect_metadata <- Transects %>%
left_join(Metadata, by = c("Site"))

#Summary of species from transects
Transect_sum <- Transect_metadata %>%
  filter(Method == "Transect") %>% #Select only data taken with visual transect method
  group_by(Site, Latitude, Census, Species) %>%
  summarise(Total_abundance = sum(Individuals))

Calculate abundances

Calculate abundances by species

First, change data frame to the wide format

Transect_wide <- Transect_sum %>%
   ungroup()  %>%
  pivot_wider(names_from = Species, values_from = Total_abundance,
              values_fill = 0)

Convert back to long format with Species and Abundance as columns. This will add rows with zeroes for all the species missing in each video

Transect_long <- Transect_wide %>%
  gather(Species, Abundance, Anisotremus_scapularis:Labrisomus_philippii)

#Obtain mean abundance per site and method
Transect_mean <- Transect_long %>%
      group_by(Site, Latitude, Species) %>%
     summarise(Mean_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))

Calculate abundances by group

Follow the same steps as abundances by type

TransectGroup_sum <- Transect_metadata %>%
  group_by(Site, Latitude, Census, Group) %>%
  summarise(Total_abundance = sum(Individuals))

TransectGroup_wide <- TransectGroup_sum %>%
  ungroup()  %>%
  pivot_wider(names_from = Group, values_from = Total_abundance,
              values_fill = 0) 
TransectGroup_wide <- TransectGroup_wide%>%
  mutate(Herbivores = Herbivore_fish + Herbivore_gastropod + Urchin + Omnivorous_crab, Invertebrate_predators = Predatory_gastropod+ Snail + Predatory_crab + Sea_star + Omnivorous_crab + Schrimp)


TransectGroup_long <- TransectGroup_wide %>%
  gather(Group, Abundance, Herbivore_fish:Invertebrate_predators)

TransectGroup_mean <- TransectGroup_long %>%
  group_by(Site, Latitude, Group) %>%
  summarise(Total_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))

Calculate Species Richness

First, add a column of presence when abundance is higher than 0.

Transect_long <- Transect_long %>%
  mutate(Presence = if_else(Abundance > 0,1,0))

Calculate sum of all species in the same site and census.

Richness_transect <- Transect_long %>%
  group_by(Site, Latitude, Census) %>%
  summarise(Richness = sum(Presence))

#Calculate mean by site  
Richness_mean <- Richness_transect %>%
  group_by(Site, Latitude)%>%
  summarise(Mean_SR = mean(Richness),n = n(), sd = sd(Richness), se = sd/sqrt(n))

Add abundances to consumption data

Calculate means and change it to the wide format

#Per Group
Group_abundances_mean <- TransectGroup_long%>%
  filter(Group == "Invertebrate_predators" | Group == "Predatory_fish"| Group == "Herbivores")%>%
  group_by(Site, Group)%>%
  summarise(Mean_Abundance = mean(Abundance)) %>%
  pivot_wider(names_from = Group, values_from = Mean_Abundance,
              values_fill = 0)

Then, add abundances data to consumption

Predation_env + group abundance summarised data

#Select only the needed columns
Group_abundances_mean_pred <- dplyr::select(Group_abundances_mean, Site, Invertebrate_predators, Predatory_fish)

Predation_4 <- Predation_3 %>%
    left_join(Group_abundances_mean_pred, by = "Site")

Herbivory_env + Benthic summarised data

#Select only the needed columns
Group_abundances_mean_herb <- dplyr::select(Group_abundances_mean, Site, Herbivores)

Herbivory_3 <- Herbivory_2 %>%
  left_join(Group_abundances_mean_herb, by = "Site")

Add species richness to consumption data

Calculate means and change it to the wide format

Richness_meanpred <- Richness_transect %>%
  group_by(Site)%>%
  summarise(Mean_SR = mean(Richness))

Then, add abundances data to consumption

Predation_4 + species richness summarised data

Predation_5 <- Predation_4 %>%
    left_join(Richness_meanpred, by = "Site")

Herbivory_3 + species richness summarised data

Herbivory_4 <- Herbivory_3 %>%
  left_join(Richness_meanpred, by = "Site")

2e): Community composition

Summarise abundances of species per site and standardise

Transect_sum <- Transect_long %>%
  group_by(Site, Latitude, Species) %>%
  summarise(Total_abundance = mean(Abundance))

## `summarise()` has grouped output by 'Site', 'Latitude'. You can override using
## the `.groups` argument.

Transect_wide2 <- Transect_sum %>%
ungroup() %>%
pivot_wider(names_from = Species, values_from = Total_abundance, values_fill = 0)

Transect_val <- Transect_wide2 %>%
  dplyr::select(- Site, -Latitude) 

Transect_stand <- decostand(Transect_val, method = "total")

Calculate Bray-Curtis dissimilariry

Transect_comp <- vegdist(Transect_stand, "bray")

Run Principal Coordinate Analysis (PcoA) to obtain the component explaining the highers variation

pcoa <- cmdscale(Transect_comp, eig = TRUE, add = TRUE)
pcoa$eig

##  [1]  1.740530e+00  1.329554e+00  8.438524e-01  6.203963e-01  5.740214e-01
##  [6]  3.517972e-01  1.671264e-01  7.777519e-02  5.878784e-02  3.878297e-02
## [11] -7.537429e-17 -2.615278e-16

(pcoa$eig/sum(pcoa$eig))*100 #proportion of variation

##  [1]  2.999557e+01  2.291298e+01  1.454260e+01  1.069165e+01  9.892446e+00
##  [6]  6.062727e+00  2.880187e+00  1.340345e+00  1.013125e+00  6.683695e-01
## [11] -1.298969e-15 -4.507061e-15

#convert pcoa results into data frame that can be plotted
pcoa_df <- data.frame(pcoa$points)
colnames(pcoa_df) <- c("PCo1", "PCo2")
pcoa_df <- pcoa_df %>%
  dplyr::select(-c(PCo2))
pcoa_df$Site <- factor(Transect_wide2$Site)
pcoa_df$Latitude <- Transect_wide2$Latitude

PCo 1 represents 30% of the variation community composition

Add PCos to consumption data

Predation_4 + PCo

pcoa_df1<- dplyr::select(pcoa_df, -Latitude)
Predation_6 <- Predation_5 %>%
    left_join(pcoa_df1, by = "Site")

Herbivory_3 + PCo

Herbivory_5 <- Herbivory_4  %>%
    left_join(pcoa_df1, by = "Site")

2f): Phlorotannins

Add metadata

Phlorotannins_metadata <- Phlorotannins %>%
left_join(Metadata, by = c("Site"))

Summarise phlorotannins per site

Phlorotannins_mean <- Phlorotannins_metadata %>%
  group_by(Site, Latitude) %>%
  summarise(phloro_mean = mean(Phlorotannins), n = n(), sd = sd(Phlorotannins), se = sd/sqrt(n))

Herbivory + phlorotannins

Phlorotannins_sum <- Phlorotannins %>%
  group_by(Site) %>%
  summarise(phloro_mean = mean(Phlorotannins))

Herbivory_6 <- Herbivory_5  %>%
    left_join(Phlorotannins_sum, by = "Site")

Part 3: Latitudinal patterns in explanatory variables

Latitudinal patterns were analysed for the recorded explanatory variables using linear models.

3a): Explore latitudinal relationships

Set general theme for plots

Theme1 <- theme_bw() + 
    theme(panel.grid.major = element_line(colour = "white"),
          panel.grid.minor = element_line(colour = "white"), axis.text.y = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 9), axis.title.x =     element_text(size = 9))

For these analyses we used the approach and code of model selection for environmental and spatial gradients provided by Anderson et al. (2022)

1) Temperature vs latitude

a. Analyses

For future plotting

Environmental_2$Latitude <- Environmental_2$Latitude*-1
#change sign of latitude

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$Temperature, xlab = "Latitude", ylab = "Temperature", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$Temperature, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$Temperature, discrete = FALSE)

## summary statistics
## ------
## min:  10.2   max:  18.7 
## median:  16 
## mean:  15.32917 
## estimated sd:  2.365349 
## estimated skewness:  -0.7647711 
## estimated kurtosis:  2.651633

normal_dist <- fitdist(Environmental_2$Temperature, "weibull")
plot(normal_dist)

T_mod1 <- glmer(Temperature ~ Latitude+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
T_mod2 <- glmer(Temperature ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
T_mod3 <- glmer(Temperature ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))

AICc(T_mod1, T_mod3, T_mod2)

##        df     AICc
## T_mod1  4 80.68285
## T_mod3  6 84.78450
## T_mod2  5 82.55903

r2(T_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.930
##      Marginal R2: 0.713

r2(T_mod2)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.935
##      Marginal R2: 0.744

r2(T_mod3)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.937
##      Marginal R2: 0.761

T_mod1 is the best fitted model

res <- resid(T_mod1)
#produce residual vs. fitted plot
plot(fitted(T_mod1), res)

#add a horizontal line at 0 
abline(0,0)

#create Q-Q plot for residuals
qqnorm(res)

#add a straight diagonal line to the plot
qqline(res)

summary(T_mod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Temperature ~ Latitude + (1 | Site)
##    Data: Environmental_2
## 
##      AIC      BIC   logLik deviance df.resid 
##     79.8     87.2    -35.9     71.8       44 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.34203 -0.41973  0.05954  0.42096  1.54302 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.004483 0.06696 
##  Residual             0.001452 0.03811 
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error t value Pr(>|z|)    
## (Intercept)  3.23945    0.20465  15.829  < 2e-16 ***
## Latitude    -0.01755    0.00669  -2.623  0.00871 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.975

r2(T_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.930
##      Marginal R2: 0.713

b. Plot predicted relationship

predictions <- ggpredict(T_mod1, terms ="Latitude")

# Plot the GAM model predictions
Temperature_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=Temperature), width = 0.08, height = 0.03, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.high, ymax = conf.low), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y =predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Temperature (°C)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

2) Salinity vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$Salinity_ppt, xlab = "Latitude", ylab = "Salinity", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$Salinity_ppt, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$Salinity_ppt, discrete = FALSE)

## summary statistics
## ------
## min:  32.8   max:  34.4 
## median:  33.8 
## mean:  33.73125 
## estimated sd:  0.3882071 
## estimated skewness:  -0.9322086 
## estimated kurtosis:  3.328625

normal_dist <- fitdist(Environmental_2$Salinity_ppt, "weibull")
plot(normal_dist)

S_mod1 <- glmer(Salinity_ppt ~ Latitude+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0215438 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

S_mod2 <- glmer(Salinity_ppt ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))
S_mod3 <- glmer(Salinity_ppt ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = Gamma(link = "log"))

check_model(S_mod1)

AICc(S_mod1, S_mod3, S_mod2)

##        df     AICc
## S_mod1  4 20.65664
## S_mod3  6 24.12187
## S_mod2  5 21.50425

BIC(S_mod1, S_mod3, S_mod2)

##        df      BIC
## S_mod1  4 27.21121
## S_mod3  6 33.30029
## S_mod2  5 29.43168

summary(S_mod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Salinity_ppt ~ Latitude + (1 | Site)
##    Data: Environmental_2
## 
##      AIC      BIC   logLik deviance df.resid 
##     19.7     27.2     -5.9     11.7       44 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8003 -0.5564  0.1690  0.5914  1.4235 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Site     (Intercept) 2.189e-05 0.004679
##  Residual             5.926e-05 0.007698
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##               Estimate Std. Error t value Pr(>|z|)    
## (Intercept)  3.5487577  0.0083718 423.896  < 2e-16 ***
## Latitude    -0.0010191  0.0002757  -3.696 0.000219 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.953
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0215438 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

r2(S_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.547
##      Marginal R2: 0.379

b. Plot predicted relationship

predictions <- ggpredict(S_mod1, terms ="Latitude [all]")

# Plot the GAM model predictions
Salinity_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=Salinity_ppt), width = 0.08, height = 0.05, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data =predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Salinity (ppt)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

3) Oxygen vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, xlab = "Latitude", ylab = "DO", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$D_oxygen_mgL, discrete = FALSE)

## summary statistics
## ------
## min:  1.62   max:  12.15 
## median:  6.445 
## mean:  6.381042 
## estimated sd:  2.626742 
## estimated skewness:  0.3329725 
## estimated kurtosis:  2.893175

normal_dist <- fitdist(Environmental_2$D_oxygen_mgL, "norm")
plot(normal_dist)

Ox_mod1 <- glmer(D_oxygen_mgL ~ Latitude+ (1|Site),  data= Environmental_2, family = gaussian)

## Warning in glmer(D_oxygen_mgL ~ Latitude + (1 | Site), data = Environmental_2,
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly

Ox_mod2 <- glmer(D_oxygen_mgL ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = gaussian)

## Warning in glmer(D_oxygen_mgL ~ poly(Latitude, 2) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly

Ox_mod3 <- glmer(D_oxygen_mgL ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = gaussian)

## Warning in glmer(D_oxygen_mgL ~ poly(Latitude, 3) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly

AICc(Ox_mod1, Ox_mod2, Ox_mod3)

##         df     AICc
## Ox_mod1  4 221.6313
## Ox_mod2  5 209.4476
## Ox_mod3  6 207.2213

BIC(Ox_mod1, Ox_mod2, Ox_mod3)

##         df      BIC
## Ox_mod1  4 228.1858
## Ox_mod2  5 217.3750
## Ox_mod3  6 216.3998

Improvement cubic model shows the best fit

Two options may be possible based on AICc. Very low df in Ox_mod 3

summary(Ox_mod3)

## Linear mixed model fit by REML ['lmerMod']
## Formula: D_oxygen_mgL ~ poly(Latitude, 3) + (1 | Site)
##    Data: Environmental_2
## 
## REML criterion at convergence: 193.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.37338 -0.33408  0.05551  0.26713  3.04460 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 4.611    2.147   
##  Residual             3.029    1.741   
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          6.3810     0.6688   9.540
## poly(Latitude, 3)1   0.6637     4.6339   0.143
## poly(Latitude, 3)2  -6.5554     4.6339  -1.415
## poly(Latitude, 3)3  -0.1882     4.6339  -0.041
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 0.000                 
## ply(Ltt,3)2 0.000  0.000          
## ply(Ltt,3)3 0.000  0.000   0.000

#Plot the GAM model predictions
Oxygen_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=D_oxygen_mgL), width = 0.08, height = 0.07, size = 1, alpha = 0.6)  + 
  Theme1 +
labs( y="Dissolved oxygen (mg"~L^-1*")", x="Latitude (°S)") +
theme(axis.title.x=element_text(size = 10),
        axis.text.x=element_text(size = 10), axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 8))

par(mfrow = c(2, 2))
check_model(Ox_mod3)

4) Visibility vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Environmental_2$Latitude, Environmental_2$V_Visibility, xlab = "Latitude", ylab = "Visibility", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Environmental_2$Latitude, Environmental_2$V_Visibility, pch = 21, col = mygrey, bg = mygrey)

descdist(Environmental_2$V_Visibility, discrete = FALSE)

## summary statistics
## ------
## min:  3   max:  10 
## median:  5.75 
## mean:  5.854167 
## estimated sd:  1.994563 
## estimated skewness:  0.4867126 
## estimated kurtosis:  2.512658

normal_dist <- fitdist(Environmental_2$V_Visibility, "norm")
plot(normal_dist)

V_mod1 <- glmer(V_Visibility ~ Latitude+ (1|Site),  data= Environmental_2, family = gaussian())

## Warning in glmer(V_Visibility ~ Latitude + (1 | Site), data = Environmental_2,
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly

V_mod2 <- glmer(V_Visibility ~poly(Latitude,2)+ (1|Site),  data= Environmental_2, family = gaussian())

## Warning in glmer(V_Visibility ~ poly(Latitude, 2) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly

V_mod3 <- glmer(V_Visibility ~poly(Latitude,3)+ (1|Site),  data= Environmental_2, family = gaussian())

## Warning in glmer(V_Visibility ~ poly(Latitude, 3) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly

check_model(V_mod1)

AICc(V_mod1, V_mod2, V_mod3)

##        df     AICc
## V_mod1  4 171.7591
## V_mod2  5 162.0855
## V_mod3  6 158.9580

BIC(V_mod1, V_mod2, V_mod3)

##        df      BIC
## V_mod1  4 178.3137
## V_mod2  5 170.0129
## V_mod3  6 168.1365

V_mod3 is the best fitted

summary(V_mod3)

## Linear mixed model fit by REML ['lmerMod']
## Formula: V_Visibility ~ poly(Latitude, 3) + (1 | Site)
##    Data: Environmental_2
## 
## REML criterion at convergence: 144.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.83174 -0.54808  0.07988  0.44641  2.45295 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 3.2862   1.8128  
##  Residual             0.8715   0.9336  
## Number of obs: 48, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          5.8542     0.5404  10.834
## poly(Latitude, 3)1   5.0601     3.7438   1.352
## poly(Latitude, 3)2  -0.1325     3.7438  -0.035
## poly(Latitude, 3)3  -4.2251     3.7438  -1.129
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 0.000                 
## ply(Ltt,3)2 0.000  0.000          
## ply(Ltt,3)3 0.000  0.000   0.000

# Plot the GAM model predictions
Visibility_plot <- ggplot() +
  geom_jitter(data = Environmental_2, aes(x= Latitude, y=V_Visibility), width = 0.08, height = 0.08, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Visibility (m)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank())

5) PCOA vs latitude

a. Analyses

pcoa_df$Latitude <- pcoa_df$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(pcoa_df$Latitude, pcoa_df$PCo1, xlab = "Latitude", ylab = "PCoA1", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(pcoa_df$Latitude, pcoa_df$PCo1, pch = 21, col = mygrey, bg = mygrey)

descdist(pcoa_df$PCo1, discrete = FALSE)

## summary statistics
## ------
## min:  -0.3828236   max:  0.6887275 
## median:  -0.1217363 
## mean:  5.247539e-17 
## estimated sd:  0.3977814 
## estimated skewness:  0.9122292 
## estimated kurtosis:  2.331104

normal_dist <- fitdist(log(pcoa_df$PCo1+1), "norm")
plot(normal_dist)

PCO_mod1 <- glm(log(PCo1 +1) ~ Latitude,  data= pcoa_df, family = gaussian())
PCO_mod2 <- glm(log(PCo1 +1) ~poly(Latitude,2),  data= pcoa_df, family = gaussian())
PCO_mod3 <- glm(log(PCo1 +1) ~poly(Latitude,3),  data= pcoa_df, family = gaussian())

check_model(PCO_mod2)

check_model(PCO_mod1)

AICc(PCO_mod1,PCO_mod2, PCO_mod3)

##          df     AICc
## PCO_mod1  3 14.72754
## PCO_mod2  4 17.28033
## PCO_mod3  5 23.21428

summary(PCO_mod1)

## 
## Call:
## glm(formula = log(PCo1 + 1) ~ Latitude, family = gaussian(), 
##     data = pcoa_df)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.85098    0.43453  -1.958   0.0786 .
## Latitude     0.02632    0.01421   1.853   0.0936 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1132384)
## 
##     Null deviance: 1.5211  on 11  degrees of freedom
## Residual deviance: 1.1324  on 10  degrees of freedom
## AIC: 11.728
## 
## Number of Fisher Scoring iterations: 2

par(mfrow = c(2, 2))
gam.check(PCO_mod1)

## 
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.

No significant effect

# Plot the GAM model predictions
PCOA_plot <- ggplot() +
  geom_jitter(data = pcoa_df, aes(x= Latitude, y=PCo1), width = 0, height = 0, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="log(PCoA1 + 1)", x="Latitude (°S)")

6) Kelp vs latitude

a. Analyses

For further plotting

Benthic_metadata_wide$Latitude <- Benthic_metadata_wide$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, xlab = "Latitude", ylab = "Kelp", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, pch = 21, col = mygrey, bg = mygrey)

Kelp_mod1 <- glmmTMB(Kelp ~ Latitude+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
Kelp_mod2 <- glmmTMB(Kelp ~poly(Latitude,2)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
Kelp_mod3 <- glmmTMB(Kelp ~poly(Latitude,3)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())

check_model(Kelp_mod3)

## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(Kelp_mod2, Kelp_mod1, Kelp_mod3)

##           df     AICc
## Kelp_mod2  6 470.8080
## Kelp_mod1  5 468.5267
## Kelp_mod3  7 473.2009

BIC(Kelp_mod2, Kelp_mod1, Kelp_mod3)

##           df      BIC
## Kelp_mod2  6 483.7652
## Kelp_mod1  5 479.4769
## Kelp_mod3  7 488.0979

summary(Kelp_mod1)

##  Family: tweedie  ( log )
## Formula:          Kelp ~ Latitude + (1 | Site)
## Data: Benthic_metadata_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    467.7    479.5   -228.8    457.7       73 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 1.159    1.076   
## Number of obs: 78, groups:  Site, 12
## 
## Dispersion parameter for tweedie family (): 6.27 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  6.42602    1.70405   3.771 0.000163 ***
## Latitude    -0.15074    0.06016  -2.506 0.012217 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2(Kelp_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.800
##      Marginal R2: 0.371

b. Plot predicted relationship

predictions <- ggpredict(Kelp_mod1, terms ="Latitude [all]")

# Plot the GAM model predictions
Kelp_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Kelp), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Kelp cover (%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

7) Bare_rock vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, xlab = "Latitude", ylab = "Bare_rock", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, pch = 21, col = mygrey, bg = mygrey)

BR_mod1 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma= 1.4, data = Benthic_metadata_wide, family = gaussian())
BR_mod2 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma=1.4, data = Benthic_metadata_wide, family = tw())


descdist(Benthic_metadata_wide$Bare_rock, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  100 
## median:  23.85013 
## mean:  36.00036 
## estimated sd:  32.50373 
## estimated skewness:  0.5844581 
## estimated kurtosis:  1.884422

normal_dist <- fitdist(Benthic_metadata_wide$Bare_rock, "exp")
plot(normal_dist)

BR_mod1 <- glmmTMB(Bare_rock ~ Latitude+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
BR_mod2 <- glmmTMB(Bare_rock ~poly(Latitude,2)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
BR_mod3 <- glmmTMB(Bare_rock ~poly(Latitude,3)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())

check_model(BR_mod2)

## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(BR_mod1, BR_mod2, BR_mod3)

##         df     AICc
## BR_mod1  5 634.5401
## BR_mod2  6 634.5044
## BR_mod3  7 636.9098

BIC(BR_mod1, BR_mod2, BR_mod3)

##         df      BIC
## BR_mod1  5 645.4903
## BR_mod2  6 647.4616
## BR_mod3  7 651.8068

summary(BR_mod2)

##  Family: tweedie  ( log )
## Formula:          Bare_rock ~ poly(Latitude, 2) + (1 | Site)
## Data: Benthic_metadata_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    633.3    647.5   -310.7    621.3       72 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.8647   0.9299  
## Number of obs: 78, groups:  Site, 12
## 
## Dispersion parameter for tweedie family ():  4.1 
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           3.126      0.283  11.046   <2e-16 ***
## poly(Latitude, 2)1    2.890      2.449   1.180    0.238    
## poly(Latitude, 2)2   -3.909      2.433  -1.606    0.108    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2(BR_mod2)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.859
##      Marginal R2: 0.225

R2 of mod2 is higher

# Plot the GAM model predictions
Bare_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Bare_rock), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Bare rock cover(%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

8) Understorey_algae vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Understorey_algae, xlab = "Latitude", ylab = "Understorey_algae", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Understorey_algae, pch = 21, col = mygrey, bg = mygrey)

descdist(Benthic_metadata_wide$Understorey_algae, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  82.92683 
## median:  0 
## mean:  10.92887 
## estimated sd:  19.35412 
## estimated skewness:  2.110129 
## estimated kurtosis:  7.04205

normal_dist <- fitdist(Benthic_metadata_wide$Understorey_algae, "exp")
plot(normal_dist)

UA_mod1 <- glmmTMB(Understorey_algae ~ Latitude+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
UA_mod2 <- glmmTMB(Understorey_algae ~poly(Latitude,2)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())
UA_mod3 <- glmmTMB(Understorey_algae ~poly(Latitude,3)+ (1|Site),  data= Benthic_metadata_wide, family = tweedie())

check_model(UA_mod2)

## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(UA_mod1, UA_mod2, UA_mod3)

##         df     AICc
## UA_mod1  5 342.5629
## UA_mod2  6 344.5077
## UA_mod3  7 342.1709

BIC(UA_mod1, UA_mod2, UA_mod3)

##         df      BIC
## UA_mod1  5 353.5131
## UA_mod2  6 357.4649
## UA_mod3  7 357.0679

summary(UA_mod1)

##  Family: tweedie  ( log )
## Formula:          Understorey_algae ~ Latitude + (1 | Site)
## Data: Benthic_metadata_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    341.7    353.5   -165.9    331.7       73 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.819    2.412   
## Number of obs: 78, groups:  Site, 12
## 
## Dispersion parameter for tweedie family (): 6.61 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  1.75831    3.33416   0.527    0.598
## Latitude    -0.04022    0.11085  -0.363    0.717

r2(UA_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.934
##      Marginal R2: 0.011

# Plot the GAM model predictions
UA_plot <- ggplot() +
  geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Understorey_algae), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="UA cover(%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

9) Predatory fish abundance vs latitude

TransectGroup_wide$Latitude <- TransectGroup_wide$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, xlab = "Latitude", ylab = "PF", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, pch = 21, col = mygrey, bg = mygrey)

descdist(TransectGroup_wide$Predatory_fish, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  333 
## median:  29 
## mean:  51.61538 
## estimated sd:  75.31432 
## estimated skewness:  2.720417 
## estimated kurtosis:  10.84418

normal_dist <- fitdist(TransectGroup_wide$Predatory_fish, "exp")
plot(normal_dist)

PF_mod1 <- glmer.nb(Predatory_fish ~ Latitude+ (1|Site),  data= TransectGroup_wide)
PF_mod2 <- glmer.nb(Predatory_fish ~poly(Latitude,2)+ (1|Site),  data= TransectGroup_wide)
PF_mod3 <- glmer.nb(Predatory_fish ~poly(Latitude,3)+ (1|Site),  data= TransectGroup_wide)

## boundary (singular) fit: see help('isSingular')

check_model(PF_mod3)

AICc(PF_mod2, PF_mod1, PF_mod3)

##         df     AICc
## PF_mod2  5 257.1651
## PF_mod1  4 257.3351
## PF_mod3  6 249.9318

BIC(PF_mod2, PF_mod1, PF_mod3)

##         df      BIC
## PF_mod2  5 260.4556
## PF_mod1  4 260.4627
## PF_mod3  6 253.0594

mod 3 is singular, so mod 1 was chosen

summary(PF_mod2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(1.2985)  ( log )
## Formula: Predatory_fish ~ poly(Latitude, 2) + (1 | Site)
##    Data: TransectGroup_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    254.2    260.5   -122.1    244.2       21 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0949 -0.6533 -0.2637  0.1770  2.9381 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4711   0.6864  
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          3.3709     0.3024  11.147  < 2e-16 ***
## poly(Latitude, 2)1  -4.8871     1.5132  -3.230  0.00124 ** 
## poly(Latitude, 2)2  -2.7915     1.5243  -1.831  0.06705 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,2)1
## ply(Ltt,2)1 0.071         
## ply(Ltt,2)2 0.159  0.175

r2(PF_mod2)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.749
##      Marginal R2: 0.546

b. Plot predicted relationship

predictions <- ggpredict(PF_mod2, terms ="Latitude [all]")

# Plot the GAM model predictions
PredFish_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Predatory_fish), width = 0.12, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Predatory Fish", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

10) Herbivorous fish abundance vs latitude

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, xlab = "Latitude", ylab = "Herbivores", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, pch = 21, col = mygrey, bg = mygrey)

descdist(TransectGroup_wide$Herbivores, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  1912 
## median:  155.5 
## mean:  389.4231 
## estimated sd:  547.3004 
## estimated skewness:  1.702098 
## estimated kurtosis:  5.106245

normal_dist <- fitdist(TransectGroup_wide$Herbivores, "exp")
plot(normal_dist)

Herb_mod1 <- glmer.nb(Herbivores ~ Latitude+ (1|Site),  data= TransectGroup_wide)
Herb_mod2 <- glmer.nb(Herbivores ~poly(Latitude,2)+ (1|Site),  data= TransectGroup_wide)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0495445 (tol = 0.002, component 1)

Herb_mod3 <- glmer.nb(Herbivores ~poly(Latitude,3)+ (1|Site),  data= TransectGroup_wide)

check_model(Herb_mod1)

AICc(Herb_mod1, Herb_mod2, Herb_mod3)

##           df     AICc
## Herb_mod1  4 326.4081
## Herb_mod2  5 326.9362
## Herb_mod3  6 327.9244

BIC(Herb_mod1, Herb_mod2, Herb_mod3)

##           df      BIC
## Herb_mod1  4 329.5358
## Herb_mod2  5 330.2266
## Herb_mod3  6 331.0519

summary(Herb_mod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(1.5227)  ( log )
## Formula: Herbivores ~ Latitude + (1 | Site)
##    Data: TransectGroup_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    324.5    329.5   -158.3    316.5       22 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2022 -0.4869 -0.0675  0.3265  1.3791 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.748    2.398   
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  10.1339     3.3117   3.060  0.00221 **
## Latitude     -0.2029     0.1101  -1.844  0.06523 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.976

r2(Herb_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.939
##      Marginal R2: 0.250

# Plot the GAM model predictions
Herb_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Herbivores), width = 0.2, height = 0.5, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Herbivores", x="Latitude (°S)")

11) Invertebrate predators abundance vs latitude

# First, examine a simple scatter plot (always a good idea!)
   plot(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, xlab = "Latitude", ylab = "Invertebrate pred", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, pch = 21, col = mygrey, bg = mygrey)

descdist(TransectGroup_wide$Invertebrate_predators, discrete = FALSE)

## summary statistics
## ------
## min:  8   max:  753 
## median:  174.5 
## mean:  242.3077 
## estimated sd:  208.9332 
## estimated skewness:  1.010871 
## estimated kurtosis:  3.417233

normal_dist <- fitdist(TransectGroup_wide$Invertebrate_predators, "exp")
plot(normal_dist)

Inv_mod1 <- glmer.nb(Invertebrate_predators ~ Latitude+ (1|Site),  data= TransectGroup_wide)
Inv_mod2 <- glmer.nb(Invertebrate_predators ~poly(Latitude,2)+ (1|Site),  data= TransectGroup_wide)
Inv_mod3 <- glmer.nb(Invertebrate_predators ~poly(Latitude,3)+ (1|Site),  data= TransectGroup_wide)

check_model(Inv_mod1)

AICc(Inv_mod2, Inv_mod1, Inv_mod3)

##          df     AICc
## Inv_mod2  5 337.9285
## Inv_mod1  4 337.7932
## Inv_mod3  6 337.1040

BIC(Inv_mod2, Inv_mod1, Inv_mod3)

##          df      BIC
## Inv_mod2  5 341.2190
## Inv_mod1  4 340.9208
## Inv_mod3  6 340.2315

summary(Inv_mod3)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Negative Binomial(3.6522)  ( log )
## Formula: Invertebrate_predators ~ poly(Latitude, 3) + (1 | Site)
##    Data: TransectGroup_wide
## 
##      AIC      BIC   logLik deviance df.resid 
##    332.7    340.2   -160.3    320.7       20 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.48939 -0.60016 -0.01322  0.64400  1.52814 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4727   0.6875  
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          5.0982     0.2268  22.483   <2e-16 ***
## poly(Latitude, 3)1  -1.1682     1.1803  -0.990   0.3223    
## poly(Latitude, 3)2  -2.3212     1.1633  -1.995   0.0460 *  
## poly(Latitude, 3)3  -2.6746     1.1575  -2.311   0.0208 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 -0.067                
## ply(Ltt,3)2  0.072 -0.069         
## ply(Ltt,3)3 -0.017  0.024  -0.085

r2(Inv_mod3)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.808
##      Marginal R2: 0.437

predictions <- ggpredict(Inv_mod3, terms ="Latitude [all]")

# Plot the GAM model predictions
Inv_plot <- ggplot() +
  geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Invertebrate_predators), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
    geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Invertebrate predators", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y = element_text(size = 7))

12) Species richness vs latitude

Richness_transect$Latitude <- Richness_transect$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Richness_transect$Latitude, Richness_transect$Richness, xlab = "Latitude", ylab = "Richness", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Richness_transect$Latitude, Richness_transect$Richness, pch = 21, col = mygrey, bg = mygrey)

descdist(Richness_transect$Richness, discrete = FALSE)

## summary statistics
## ------
## min:  2   max:  15 
## median:  8 
## mean:  8.115385 
## estimated sd:  3.076712 
## estimated skewness:  0.3378056 
## estimated kurtosis:  3.339398

normal_dist <- fitdist(Richness_transect$Richness, "norm")
plot(normal_dist)

R_mod1 <- glmer(Richness ~ Latitude+ (1|Site),  data= Richness_transect, family = gaussian())

## Warning in glmer(Richness ~ Latitude + (1 | Site), data = Richness_transect, :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly

R_mod2 <- glmer(Richness ~poly(Latitude,2)+ (1|Site),  data= Richness_transect, family = gaussian())

## Warning in glmer(Richness ~ poly(Latitude, 2) + (1 | Site), data =
## Richness_transect, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly

R_mod3 <- glmer(Richness ~poly(Latitude,3)+ (1|Site),  data= Richness_transect, family = gaussian())

## Warning in glmer(Richness ~ poly(Latitude, 3) + (1 | Site), data =
## Richness_transect, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly

AICc(R_mod2, R_mod1, R_mod3)

##        df     AICc
## R_mod2  5 127.3899
## R_mod1  4 136.8463
## R_mod3  6 124.1031

BIC(R_mod2, R_mod1, R_mod3)

##        df      BIC
## R_mod2  5 130.6804
## R_mod1  4 139.9739
## R_mod3  6 127.2306

summary(R_mod3)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Richness ~ poly(Latitude, 3) + (1 | Site)
##    Data: Richness_transect
## 
## REML criterion at convergence: 107.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6252 -0.5226  0.0082  0.3615  2.0845 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 5.464    2.338   
##  Residual             4.171    2.042   
## Number of obs: 26, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          8.0744     0.7901  10.220
## poly(Latitude, 3)1  -3.4362     4.0846  -0.841
## poly(Latitude, 3)2  -3.0754     4.0715  -0.755
## poly(Latitude, 3)3  -5.9172     4.0749  -1.452
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 -0.079                
## ply(Ltt,3)2  0.057 -0.081         
## ply(Ltt,3)3 -0.014  0.039  -0.072

r2(R_mod3)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.649
##      Marginal R2: 0.189

No significant effect

# Plot the GAM model predictions
Richness_plot <- ggplot() +
  geom_jitter(data = Richness_transect, aes(x= Latitude, y=Richness), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  Theme1 +
labs( y="Species Richness", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y = element_text(size = 8))

13) Fish length vs latitude

Size_fish <- filter(Size_mean_sp, Type == "Fish")
Size_fish$Latitude <- Size_fish$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Size_fish$Latitude, Size_fish$mean_Size_sp, xlab = "Latitude", ylab = "Fish Length", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Size_fish$Latitude, Size_fish$mean_Size_sp, pch = 21, col = mygrey, bg = mygrey)

descdist(Size_fish$mean_Size_sp, discrete = FALSE)

## summary statistics
## ------
## min:  5.4   max:  65 
## median:  23 
## mean:  22.49199 
## estimated sd:  9.684538 
## estimated skewness:  0.9708081 
## estimated kurtosis:  5.717812

normal_dist <- fitdist(Size_fish$mean_Size_sp, "logis")
plot(normal_dist)

Length_mod1 <- glmer(mean_Size_sp ~ Latitude+ (1|Site) + (1|Species),  data= Size_fish, family = Gamma(link = "log"))
Length_mod2 <- glmer(mean_Size_sp ~poly(Latitude,2)+ (1|Site) + (1|Species),  data= Size_fish, family = Gamma(link = "log"))
Length_mod3 <- glmer(mean_Size_sp ~poly(Latitude,3)+ (1|Site) + (1|Species),  data= Size_fish, family = Gamma(link = "log"))

check_model(Length_mod1)

AICc(Length_mod1, Length_mod2, Length_mod3)

##             df     AICc
## Length_mod1  5 773.3847
## Length_mod2  6 772.3818
## Length_mod3  7 774.6382

BIC(Length_mod1, Length_mod2, Length_mod3)

##             df      BIC
## Length_mod1  5 786.7493
## Length_mod2  6 788.3065
## Length_mod3  7 793.0830

summary(Length_mod2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean_Size_sp ~ poly(Latitude, 2) + (1 | Site) + (1 | Species)
##    Data: Size_fish
## 
##      AIC      BIC   logLik deviance df.resid 
##    771.6    788.3   -379.8    759.6      113 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2255 -0.4292  0.0110  0.4239  4.0786 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Species  (Intercept) 0.091860 0.30308 
##  Site     (Intercept) 0.002929 0.05412 
##  Residual             0.073020 0.27022 
## Number of obs: 119, groups:  Species, 28; Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value Pr(>|z|)    
## (Intercept)          3.0014     0.1174  25.566   <2e-16 ***
## poly(Latitude, 2)1   0.1303     0.3940   0.331   0.7409    
## poly(Latitude, 2)2  -0.7157     0.3493  -2.049   0.0405 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,2)1
## ply(Ltt,2)1 -0.085        
## ply(Ltt,2)2 -0.027  0.036

r2(Length_mod2)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.585
##      Marginal R2: 0.026

predictions <- ggpredict(Length_mod2, terms ="Latitude [all]")

# Plot the GAM model predictions
Size_plot <- ggplot() +
  geom_jitter(data = Size_fish, aes(x= Latitude, y=mean_Size_sp), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Fish length (cm)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(), axis.title.y=element_text(8))

14) Phlorotannins vs latitude

Phlorotannins_metadata$Latitude <- Phlorotannins_metadata$Latitude*-1

# First, examine a simple scatter plot (always a good idea!)
   plot(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, xlab = "Latitude", ylab = "phlorotannins", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, pch = 21, col = mygrey, bg = mygrey)

descdist(Phlorotannins_metadata$Phlorotannins, discrete = FALSE)

## summary statistics
## ------
## min:  1.83   max:  94.69 
## median:  11.625 
## mean:  17.26067 
## estimated sd:  15.2777 
## estimated skewness:  2.275762 
## estimated kurtosis:  9.890266

normal_dist <- fitdist(Phlorotannins_metadata$Phlorotannins, "exp")
plot(normal_dist)

Phloro_mod1 <- glmer(Phlorotannins ~ Latitude+ (1|Site),  data= Phlorotannins_metadata, family = Gamma(link= "log"))
Phloro_mod2 <- glmer(Phlorotannins ~poly(Latitude,2)+ (1|Site),  data= Phlorotannins_metadata, family = Gamma(link= "log"))
Phloro_mod3 <- glmer(Phlorotannins ~poly(Latitude,3)+ (1|Site),  data= Phlorotannins_metadata, family = Gamma(link= "log"))

check_model(Phloro_mod1)

AICc(Phloro_mod1, Phloro_mod2, Phloro_mod3)

##             df     AICc
## Phloro_mod1  4 831.2629
## Phloro_mod2  5 832.2302
## Phloro_mod3  6 828.0473

summary(Phloro_mod3)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Phlorotannins ~ poly(Latitude, 3) + (1 | Site)
##    Data: Phlorotannins_metadata
## 
##      AIC      BIC   logLik deviance df.resid 
##    827.3    844.0   -407.7    815.3      114 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5999 -0.5788 -0.0462  0.2840  3.8222 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Site     (Intercept) 0.06263  0.2503  
##  Residual             0.29291  0.5412  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                    Estimate Std. Error t value Pr(>|z|)    
## (Intercept)          2.6600     0.1165  22.838  < 2e-16 ***
## poly(Latitude, 3)1  -3.1963     1.2765  -2.504  0.01228 *  
## poly(Latitude, 3)2  -2.0276     1.2753  -1.590  0.11185    
## poly(Latitude, 3)3   3.8707     1.2767   3.032  0.00243 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1  0.000                
## ply(Ltt,3)2  0.001  0.000         
## ply(Ltt,3)3  0.000 -0.001  -0.013

r2(Phloro_mod3)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.546
##      Marginal R2: 0.435

predictions <- ggpredict(Phloro_mod3, terms ="Latitude [all]")

# Plot the GAM model predictions
Phloro_plot <- ggplot() +
  geom_jitter(data = Phlorotannins_metadata, aes(x= Latitude, y=Phlorotannins), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
  geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
  geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) + 
  Theme1 +
labs( y="Phlorotannins (mg"~g^-1*" DW)", x="Latitude (°S)")

Final plots

Figure 2 - Abiotic predictors

Temperature_plot + Salinity_plot + Oxygen_plot + Phloro_plot+ plot_layout(ncol=2)

Figure 2: Latitudinal variation of abiotic variables: temperature (°C), salinity (ppt) and dissolved oxygen (mgL) and visibility (m).

Save the plot in tiff format

ggsave("Latitudevsabiotic.jpeg", units="mm", width=125, height=85, dpi=300)

Figure 3

Visibility_plot + Bare_plot+ Kelp_plot + UA_plot + Inv_plot + PredFish_plot + Size_plot + Richness_plot + PCOA_plot+ Herb_plot + plot_layout(ncol=2)

## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database

Figure 3: Latitudinal variation of biotic and habitat variables

ggsave("FigurelatvsBio.jpeg", units="mm", width=110, height=140, dpi=300)

Part 4: Latitudinal patterns in consumption

Latitudinal patterns in consumption (predation and herbivory) were explored and analysed

4a): Predation vs latitude

First change latitude to positive numbers for plotting

Predation_env$Latitude <- Predation_env$Latitude*-1

C_24 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 24)
C_2 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 2)
S_24 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 24)
S_2 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 2)

1) Crabs 2h

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(C_2$Latitude, C_2$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(C_2$Latitude, C_2$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)

descdist(C_2$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  1 
## mean:  1.958333 
## estimated sd:  2.043213 
## estimated skewness:  0.4960334 
## estimated kurtosis:  1.571122

normal_dist <- fitdist(C_2$Consumed_individuals, "exp")
plot(normal_dist)

C2_mod1 <- glmer(Consumed_individuals ~ Latitude+ (1|Site),  data= C_2, family = poisson())
C2_mod2 <- glmer(Consumed_individuals ~poly(Latitude,2)+ (1|Site),  data= C_2, family = poisson())
C2_mod3 <- glmer(Consumed_individuals ~poly(Latitude,3)+ (1|Site),  data= C_2, family = poisson())

check_model(C2_mod1)

AICc(C2_mod2, C2_mod1, C2_mod3)

##         df     AICc
## C2_mod2  4 424.9806
## C2_mod1  3 423.1983
## C2_mod3  5 425.7216

BIC(C2_mod2, C2_mod1, C2_mod3)

##         df      BIC
## C2_mod2  4 435.7828
## C2_mod1  3 431.3539
## C2_mod3  5 439.1328

summary(C2_mod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Latitude + (1 | Site)
##    Data: C_2
## 
##      AIC      BIC   logLik deviance df.resid 
##    423.0    431.4   -208.5    417.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92403 -0.94180  0.03414  0.62864  3.02889 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.414    0.6434  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  2.33581    0.91197   2.561   0.0104 *
## Latitude    -0.06553    0.03032  -2.161   0.0307 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.975

r2(C2_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.601
##      Marginal R2: 0.198

2) Crabs 24h

a. Analyses

# First, examine a simple scatter plot (always a good idea!)
   plot(C_24$Latitude, C_24$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed 24 h", las=1)
   mygrey <- grey(level = 0.65, alpha = 0.4)
   points(C_24$Latitude, C_24$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)

descdist(C_24$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  5 
## mean:  3.791667 
## estimated sd:  1.511279 
## estimated skewness:  -0.8418895 
## estimated kurtosis:  2.211133

normal_dist <- fitdist(C_24$Consumed_individuals, "exp")
plot(normal_dist)

C24_mod1 <- glm(Consumed_individuals ~ Latitude,  data= C_24, family = poisson())
C24_mod2 <- glm(Consumed_individuals ~poly(Latitude,2),  data= C_24, family = poisson())
C24_mod3 <- glm(Consumed_individuals ~poly(Latitude,3),  data= C_24, family = poisson())
check_overdispersion(C24_mod1)

## # Overdispersion test
## 
##        dispersion ratio =  0.420
##   Pearson's Chi-Squared = 49.526
##                 p-value =      1

## No overdispersion detected.

check_model(C24_mod1)

AICc(C24_mod1, C24_mod2, C24_mod3)

##          df     AICc
## C24_mod1  2 433.7401
## C24_mod2  3 434.2198
## C24_mod3  4 436.0433

BIC(C24_mod1, C24_mod2, C24_mod3)

##          df      BIC
## C24_mod1  2 439.2125
## C24_mod2  3 442.3754
## C24_mod3  4 446.8454

summary(C24_mod1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Latitude, family = poisson(), 
##     data = C_24)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.452269   0.208192  11.779  < 2e-16 ***
## Latitude    -0.038689   0.007222  -5.357 8.46e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 88.137  on 119  degrees of freedom
## Residual deviance: 58.265  on 118  degrees of freedom
## AIC: 433.64
## 
## Number of Fisher Scoring iterations: 4

r2(C24_mod1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.424

Calculate predictions for Crab consumption after 2 and 24 h

predictions <- ggpredict(C2_mod1, terms ="Latitude")
predictions1 <- ggpredict(C24_mod1, terms ="Latitude")

Add new column to mix datasets

predictions$Assay_time <- "2 h"
predictions1$Assay_time <- "24 h"
C_2$Assay_time <- "2 h"
C_24$Assay_time <- "24 h"

#Bind datasets

total <- rbind(predictions, predictions1)

total2 <- rbind (C_2, C_24)

3) Crabs 2 and 24 h prediction plot

ggplot()+
  geom_ribbon(data = total, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, group = Assay_time, fill = Assay_time ), alpha = .25)  +
  geom_line(data = total, aes(x = x, y = predicted, color = Assay_time, group = Assay_time), size = 1) +
  geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = 0.5)+
  labs(x = "Latitude", y = "Crabs consumed")  + Theme1 +
  scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
   scale_fill_manual("Assay time", values=c("darkslateblue", "orange1")) + 
  theme(legend.position = c(0.86, 0.81),legend.title = element_text(size = 8), 
               legend.text = element_text(size = 6.3), legend.background = element_rect(fill="NA",
                                         size=1, 
                                         color ="NA")) +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).

Save plot in tiff format

ggsave("Crabs.jpeg", units="mm", width=100, height=65, dpi=300)

## Warning: Removed 19 rows containing missing values or values outside the scale range
## (`geom_point()`).

4) Squidpops 2h

a. Analyses

descdist(S_2$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  2 
## mean:  2.375 
## estimated sd:  1.974895 
## estimated skewness:  0.1159755 
## estimated kurtosis:  1.430103

normal_dist <- fitdist(S_2$Consumed_individuals, "nbinom")
plot(normal_dist)

S2_mod1 <- glmer(Consumed_individuals ~ Latitude + (1|Site),  data= S_2, family = poisson())
S2_mod2 <- glmer(Consumed_individuals ~poly(Latitude,2)+ (1|Site),  data= S_2, family = poisson())
S2_mod3 <- glmer(Consumed_individuals ~poly(Latitude,3)+ (1|Site),  data= S_2, family = poisson())
check_overdispersion(S2_mod1)

## # Overdispersion test
## 
##        dispersion ratio =   1.050
##   Pearson's Chi-Squared = 122.805
##                 p-value =   0.338

## No overdispersion detected.

check_model(S2_mod1)

AICc(S2_mod2, S2_mod1, S2_mod3)

##         df     AICc
## S2_mod2  4 462.0740
## S2_mod1  3 460.3905
## S2_mod3  5 463.5723

BIC(S2_mod2, S2_mod1, S2_mod3)

##         df      BIC
## S2_mod2  4 472.8761
## S2_mod1  3 468.5461
## S2_mod3  5 476.9835

summary(S2_mod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Latitude + (1 | Site)
##    Data: S_2
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.2    468.5   -227.1    454.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7204 -0.8067 -0.0802  0.7098  2.3659 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5241   0.7239  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.14773    0.99008   0.149    0.881
## Latitude     0.01733    0.03209   0.540    0.589
## 
## Correlation of Fixed Effects:
##          (Intr)
## Latitude -0.974

r2(S2_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.618
##      Marginal R2: 0.016

5) Squidpops 24h

a. Analyses

descdist(S_24$Consumed_individuals, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  5 
## median:  5 
## mean:  4.25 
## estimated sd:  1.258584 
## estimated skewness:  -1.822846 
## estimated kurtosis:  5.495426

normal_dist <- fitdist(S_24$Consumed_individuals, "pois")
plot(normal_dist)

S24_mod1 <- glm(Consumed_individuals ~ Latitude,  data= S_24, family = poisson())
S24_mod2 <- glm(Consumed_individuals ~poly(Latitude,2),  data= S_24, family = poisson())
S24_mod3 <- glm(Consumed_individuals ~poly(Latitude,3),  data= S_24, family = poisson())
check_overdispersion(S24_mod1)

## # Overdispersion test
## 
##        dispersion ratio =  0.348
##   Pearson's Chi-Squared = 41.044
##                 p-value =      1

## No overdispersion detected.

check_model(S24_mod1)

AICc(S24_mod2, S24_mod1, S24_mod3)

##          df     AICc
## S24_mod2  3 447.7243
## S24_mod1  2 446.1367
## S24_mod3  4 449.8652

BIC(S24_mod2, S24_mod1, S24_mod3)

##          df      BIC
## S24_mod2  3 455.8799
## S24_mod1  2 451.6091
## S24_mod3  4 460.6673

summary(S24_mod1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Latitude, family = poisson(), 
##     data = S_24)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.972899   0.196634  10.033  < 2e-16 ***
## Latitude    -0.017890   0.006608  -2.707  0.00678 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 54.280  on 118  degrees of freedom
## AIC: 446.03
## 
## Number of Fisher Scoring iterations: 4

r2(S24_mod1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.149

#For squidpops 24 h 
predictions3 <- ggpredict(S24_mod1, terms ="Latitude")

Add new column to mix datasets

S_2$Assay_time <- "2 h"
S_24$Assay_time <- "24 h"
predictions3$Assay_time <- "24 h"
##Keep only colums for plotting
total2 <- rbind (S_2, S_24)

3) Squidpops 2 and 24 h prediction plot

ggplot() +
  geom_ribbon(data = predictions3, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, group = Assay_time, fill = Assay_time ), alpha = .25) +
  geom_line(data = predictions3, aes(x = x, y = predicted, color = Assay_time, group = Assay_time ), size = 1) +
  geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = .5)+
  labs(x = "Latitude", y = "Squidpops consumed")  + Theme1 +
  scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
   scale_fill_manual("Assay time", values=c("orange1")) + 
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))

## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

Save plot in tiff format

ggsave("Squidpops.jpeg", units="mm", width=100, height=65, dpi=300)

## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).

4b): Herbivory vs latitude

Herbivory_env$Latitude <- Herbivory_env$Latitude*-1

a. Analyses

descdist(Herbivory_env$Consumed_proportion, discrete = FALSE)

## summary statistics
## ------
## min:  0   max:  1 
## median:  0 
## mean:  0.2260169 
## estimated sd:  0.3698909 
## estimated skewness:  1.239303 
## estimated kurtosis:  2.780775

normal_dist <- fitdist(Herbivory_env$Consumed_proportion, "exp")
plot(normal_dist)

Herb_mod1 <- glmmTMB(Consumed_proportion ~ Latitude+ (1|Site),  data= Herbivory_env, family = tweedie())
Herb_mod2 <- glmmTMB(Consumed_proportion ~poly(Latitude,2)+ (1|Site),  data= Herbivory_env, family = tweedie())
Herb_mod3 <- glmmTMB(Consumed_proportion ~poly(Latitude,3)+ (1|Site),  data= Herbivory_env, family = tweedie())

check_model(Herb_mod1)

## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(Herb_mod2, Herb_mod1, Herb_mod3)

##           df     AICc
## Herb_mod2  6 732.3190
## Herb_mod1  5 730.7644
## Herb_mod3  7 732.0831

summary(Herb_mod1)

##  Family: tweedie  ( log )
## Formula:          Consumed_proportion ~ Latitude + (1 | Site)
## Data: Herbivory_env
## 
##      AIC      BIC   logLik deviance df.resid 
##    730.7    752.6   -360.3    720.7      585 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.69     2.586   
## Number of obs: 590, groups:  Site, 12
## 
## Dispersion parameter for tweedie family (): 0.875 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -8.9565     3.5611  -2.515   0.0119 *
## Latitude      0.1899     0.1138   1.668   0.0953 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2(Herb_mod1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.956
##      Marginal R2: 0.195

b. Plot predicted relationship

ggplot() +
  geom_jitter(data = Herbivory_env, aes(x= Latitude, y=Consumed_proportion), color = "orange1", width=0.2, height = 0.02, size = 1, alpha = .5)+
  labs(x = "Latitude (°S)", y = "Blade consumed proportion")  + Theme1  + 
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(15, 40, 5)) +
  scale_y_continuous(limits = c(-0.1,1.2), breaks = seq(0, 1, 0.2))

Save plot to tiff format

ggsave("Seaweed24.jpeg", units="mm", width=105, height=65, dpi=300)

Part 5: Predation intensity: Effect of site-specific variables

Here, the individual effects of each predictor will be tested.

First, create datasets to be used in the analysis

Predation_analysis <- Predation_6 %>%
  dplyr::select(-Rope, -Bait, -Camera, -Assay_time, -Initial_Individuals, -Remaining_individuals, -Survival_proportion,-Consumed_proportion)

#Separate by Bait and Assay time 

Crabs_24 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 24)
Crabs_2 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 2)
Squidpops_24 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 24)
Squidpops_2 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 2)

5a): Standardisation and test autocorrelation

1) Standardise

#Crabs 24h
Crabs24_num <- Crabs_24 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude, -Day_length) 
Crabs24_stand <- decostand(Crabs24_num, method="range", na.rm = TRUE)
#Test autocorrelation 
corrplot(cor(Crabs24_stand), method = "number")

No multicolinearity found

Crabs24_stand$Site <- Crabs_24$Site
Crabs24_stand$Consumed_individuals<- Crabs_24$Consumed_individuals
Crabs24_stand$Latitude<- Crabs_24$Latitude

#Crabs 2h 
Crabs2_num <- Crabs_2 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Crabs2_stand <- decostand(Crabs2_num, method="range", na.rm = TRUE)
Crabs2_stand$Site <- Crabs_2$Site
Crabs2_stand$Consumed_individuals<- Crabs_2$Consumed_individuals
Crabs2_stand$Latitude<- Crabs_2$Latitude
#Squidpops 24h
Squidpops24_num <- Squidpops_24 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Squidpops24_stand <- decostand(Squidpops24_num, method="range", na.rm = TRUE)
Squidpops24_stand$Site <- Squidpops_24$Site
Squidpops24_stand$Consumed_individuals<- Squidpops_24$Consumed_individuals
Squidpops24_stand$Latitude<- Squidpops_24$Latitude
#Squidpops 2h
Squidpops2_num <- Squidpops_2 %>%
  dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude) 
Squidpops2_stand <- decostand(Squidpops2_num, method="range", na.rm = TRUE)
Squidpops2_stand$Site <- Squidpops_2$Site
Squidpops2_stand$Consumed_individuals<- Squidpops_2$Consumed_individuals
Squidpops2_stand$Latitude<- Squidpops_2$Latitude

5b) Model selection: Crabs 2h

Temperature

my_glmT1<- glmer(Consumed_individuals~  Temperature + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmT1)

AICc(my_glmT1)

## [1] 420.2918

summary(my_glmT1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    420.1    428.4   -207.0    414.1      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93381 -0.92312  0.06015  0.65176  2.79744 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2979   0.5458  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -0.8086     0.4487  -1.802  0.07151 . 
## Temperature   1.8546     0.6104   3.038  0.00238 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.918

r2(my_glmT1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.599
##      Marginal R2: 0.307

Salinity

my_glmSC2_1<- glmer(Consumed_individuals~  Salinity_ppt + (1|Site), data=Crabs2_stand, family = poisson)


check_overdispersion(my_glmSC2_1)

## # Overdispersion test
## 
##        dispersion ratio =   1.189
##   Pearson's Chi-Squared = 139.141
##                 p-value =    0.08

## No overdispersion detected.

check_model(my_glmSC2_1)

AICc(my_glmSC2_1)

## [1] 424.5386

summary(my_glmSC2_1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    424.3    432.7   -209.2    418.3      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92849 -0.91521  0.04145  0.66421  2.99635 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.468    0.6841  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -0.6810     0.6667  -1.022   0.3070  
## Salinity_ppt   1.5821     0.9192   1.721   0.0852 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.947

r2(my_glmSC2_1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.603
##      Marginal R2: 0.148

Dissolved oxygen

my_glmOxC2_1<- glmer(Consumed_individuals~  D_oxygen_mgL + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmOxC2_1)

AICc(my_glmOxC2_1)

## [1] 422.5227

summary(my_glmOxC2_1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    422.3    430.7   -208.2    416.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9665 -0.9020  0.0611  0.5761  3.2095 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.365    0.6042  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -0.2934     0.3526  -0.832   0.4054  
## D_oxygen_mgL   1.5663     0.6451   2.428   0.0152 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.836

r2(my_glmOxC2_1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.585
##      Marginal R2: 0.215

Visibility

my_glmVC2_1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmVC2_1)

AICc(my_glmVC2_1)

## [1] 427.2626

summary(my_glmVC2_1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    427.1    435.4   -210.5    421.1      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93460 -0.89601  0.07717  0.53694  3.02786 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5924   0.7697  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.39200    0.44714   0.877    0.381
## V_Visibility -0.01098    0.89880  -0.012    0.990
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.847

r2(my_glmVC2_1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.592
##      Marginal R2: 0.000

Community composition

my_glmC2_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_PCO1)

AICc(my_glmC2_PCO1)

## [1] 426.871

summary(my_glmC2_PCO1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.7    435.0   -210.3    420.7      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92572 -0.90531  0.06451  0.54796  2.99751 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5818   0.7628  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.5359     0.3320   1.614    0.107
## PCo1         -0.4240     0.6764  -0.627    0.531
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.703

r2(my_glmC2_PCO1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.597
##      Marginal R2: 0.023

Kelp Cover

my_glmC2_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_K1)

AICc(my_glmC2_K1)

## [1] 424.3606

summary(my_glmC2_K1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    424.2    432.5   -209.1    418.2      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92113 -0.89423  0.02569  0.55099  3.16011 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4494   0.6704  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.02724    0.29757   0.092   0.9271  
## Kelp         1.32373    0.73207   1.808   0.0706 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Kelp -0.705

r2(my_glmC2_K1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.589
##      Marginal R2: 0.138

Bare rock cover

my_glmC2_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_BR1)

AICc(my_glmC2_BR1)

## [1] 424.1326

summary(my_glmC2_BR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    423.9    432.3   -209.0    417.9      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9360 -0.9040  0.0409  0.5802  3.2149 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4321   0.6573  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.8447     0.3094   2.730  0.00633 **
## Bare_rock    -1.1994     0.6315  -1.899  0.05751 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.741

r2(my_glmC2_BR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.593
##      Marginal R2: 0.163

Understorey algae cover

my_glmC2_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_UA1)

AICc(my_glmC2_UA1)

## [1] 427.0992

summary(my_glmC2_UA1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.9    435.3   -210.4    420.9      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92743 -0.90039  0.08966  0.54638  3.04101 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5829   0.7635  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)
## (Intercept)         0.4574     0.2908   1.573    0.116
## Understorey_algae  -0.2698     0.6637  -0.406    0.684
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.584

r2(my_glmC2_UA1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.592
##      Marginal R2: 0.009

Species richness

my_glmC2_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_SR1)

AICc(my_glmC2_SR1)

## [1] 425.1803

summary(my_glmC2_SR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    425.0    433.3   -209.5    419.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.90607 -0.90370  0.06207  0.59243  3.08960 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4998   0.707   
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.1651     0.4393  -0.376    0.707
## Mean_SR       1.1522     0.7741   1.488    0.137
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.864

r2(my_glmC2_SR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.598
##      Marginal R2: 0.105

Invertebrate predators

my_glmC2_Inv1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_Inv1)

AICc(my_glmC2_Inv1)

## [1] 427.1083

summary(my_glmC2_Inv1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.9    435.3   -210.5    420.9      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93478 -0.89816  0.08557  0.54475  3.04770 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5865   0.7658  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)
## (Intercept)              0.2603     0.4017   0.648    0.517
## Invertebrate_predators   0.2879     0.7302   0.394    0.693
## 
## Correlation of Fixed Effects:
##             (Intr)
## Invrtbrt_pr -0.808

r2(my_glmC2_Inv1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.593
##      Marginal R2: 0.009

Predatory fish abundance

my_glmC2_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_PF1)

AICc(my_glmC2_PF1)

## [1] 426.858

summary(my_glmC2_PF1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.7    435.0   -210.3    420.7      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92759 -0.89488  0.07857  0.55013  3.05066 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5723   0.7565  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)      0.2205     0.3528   0.625    0.532
## Predatory_fish   0.4397     0.6854   0.642    0.521
## 
## Correlation of Fixed Effects:
##             (Intr)
## Predtry_fsh -0.748

r2(my_glmC2_PF1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.593
##      Marginal R2: 0.022

Fish size

my_glmC2_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Crabs2_stand, family = poisson)

check_model(my_glmC2_FL1)

AICc(my_glmC2_FL1)

## [1] 426.9261

summary(my_glmC2_FL1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    426.7    435.1   -210.4    420.7      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93714 -0.91055  0.06116  0.54889  3.03983 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5781   0.7603  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.06379    0.60517   0.105    0.916
## Fish_size    0.52233    0.89451   0.584    0.559
## 
## Correlation of Fixed Effects:
##           (Intr)
## Fish_size -0.921

r2(my_glmC2_FL1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.594
##      Marginal R2: 0.020

5c) Model selection: Crabs 24h

Temperature

my_glmC24_T1<- glmer(Consumed_individuals~Temperature + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_T1)

AICc(my_glmC24_T1)

## [1] 438.9369

summary(my_glmC24_T1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    438.7    447.1   -216.4    432.7      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1336 -0.2824  0.1619  0.4568  1.3189 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.
##  Site   (Intercept) 0.0003175 0.01782 
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.7379     0.1353   5.456 4.88e-08 ***
## Temperature   0.8747     0.1776   4.925 8.42e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.936

r2(my_glmC24_T1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.230
##      Marginal R2: 0.229

Salinity

my_glmC24_S1<- glmer(Consumed_individuals~Salinity_ppt + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_S1)

AICc(my_glmC24_S1)

## [1] 443.1697

summary(my_glmC24_S1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    443.0    451.3   -218.5    437.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.01001 -0.46529  0.09323  0.42898  1.47247 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.01181  0.1087  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.6358     0.1945   3.269 0.001078 ** 
## Salinity_ppt   0.9897     0.2621   3.777 0.000159 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.956

r2(my_glmC24_S1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.233
##      Marginal R2: 0.195

Dissolved oxygen

my_glmC24_Ox1<- glmer(Consumed_individuals~D_oxygen_mgL + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_Ox1)

AICc(my_glmC24_Ox1)

## [1] 452.5543

summary(my_glmC24_Ox1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.3    460.7   -223.2    446.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0829 -0.3910  0.1934  0.3175  1.2599 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05924  0.2434  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.1962     0.1534   7.796 6.38e-15 ***
## D_oxygen_mgL   0.2413     0.2866   0.842      0.4    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.831

r2(my_glmC24_Ox1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.017

Visibility

my_glmC24_V1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_V1)

AICc(my_glmC24_V1)

## [1] 452.5114

summary(my_glmC24_V1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.3    460.7   -223.2    446.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0495 -0.3868  0.2021  0.3902  1.2724 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05948  0.2439  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.4210     0.1604   8.858   <2e-16 ***
## V_Visibility  -0.2825     0.3267  -0.865    0.387    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.846

r2(my_glmC24_V1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.218
##      Marginal R2: 0.019

Consumer community

my_glmC24_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_PCO1)

AICc(my_glmC24_PCO1)

## [1] 450.5558

summary(my_glmC24_PCO1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    450.3    458.7   -222.2    444.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9758 -0.3167  0.2087  0.3774  1.2274 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.04741  0.2177  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.4411     0.1109  12.993   <2e-16 ***
## PCo1         -0.3920     0.2300  -1.705   0.0883 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.698

r2(my_glmC24_PCO1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.223
##      Marginal R2: 0.065

Kelp cover

my_glmC24_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_K1)

AICc(my_glmC24_K1)

## [1] 448.7633

summary(my_glmC24_K1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    448.6    456.9   -221.3    442.6      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0196 -0.4200  0.1414  0.4214  1.2459 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.03546  0.1883  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1427     0.1040  10.987   <2e-16 ***
## Kelp          0.5844     0.2518   2.321   0.0203 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## Kelp -0.715

r2(my_glmC24_K1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.210
##      Marginal R2: 0.090

Bare rock cover

my_glmC24_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_BR1)

AICc(my_glmC24_BR1)

## [1] 447.5098

summary(my_glmC24_BR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    447.3    455.7   -220.7    441.3      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0745 -0.3637  0.1577  0.3357  1.5700 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.02761  0.1662  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.51816    0.09993  15.192  < 2e-16 ***
## Bare_rock   -0.56632    0.20938  -2.705  0.00683 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.733

r2(my_glmC24_BR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.123

Understorey algae cover

my_glmC24_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_UA1)

AICc(my_glmC24_UA1)

## [1] 453.1994

summary(my_glmC24_UA1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    453.0    461.4   -223.5    447.0      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0561 -0.4182  0.1760  0.3756  1.3128 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.06408  0.2531  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.28831    0.10876  11.846   <2e-16 ***
## Understorey_algae  0.05185    0.24524   0.211    0.833    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.590

r2(my_glmC24_UA1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.001

Species richness

my_glmC24_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_SR1)

AICc(my_glmC24_SR1)

## [1] 451.6355

summary(my_glmC24_SR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.4    459.8   -222.7    445.4      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9950 -0.3586  0.1690  0.4479  1.2775 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.05277  0.2297  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1192     0.1647   6.797 1.07e-11 ***
## Mean_SR       0.3830     0.2925   1.309     0.19    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.867

r2(my_glmC24_SR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.040

Abundance of invertebrate predators

my_glmC24_Inv1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_Inv1)

AICc(my_glmC24_Inv1)

## [1] 453.0005

summary(my_glmC24_Inv1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.8    461.2   -223.4    446.8      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0448 -0.4524  0.1776  0.4005  1.3483 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.06186  0.2487  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              1.3610     0.1462   9.308   <2e-16 ***
## Invertebrate_predators  -0.1345     0.2703  -0.498    0.619    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Invrtbrt_pr -0.805

r2(my_glmC24_Inv1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.214
##      Marginal R2: 0.007

Abundance of fish predators

my_glmC24_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_PF1)

AICc(my_glmC24_PF1)

## [1] 452.9183

summary(my_glmC24_PF1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.7    461.1   -223.4    446.7      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0344 -0.3964  0.1898  0.4233  1.2267 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.06197  0.2489  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.2463     0.1308   9.528   <2e-16 ***
## Predatory_fish   0.1465     0.2549   0.575    0.565    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Predtry_fsh -0.748

r2(my_glmC24_PF1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.216
##      Marginal R2: 0.008

Fish legth (cm)

my_glmC24_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Crabs24_stand, family = poisson)

check_model(my_glmC24_FL1)

AICc(my_glmC24_FL1)

## [1] 448.4084

summary(my_glmC24_FL1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    448.2    456.6   -221.1    442.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0671 -0.2882  0.2492  0.3076  1.0942 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.03566  0.1888  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.8921     0.1932   4.618 3.88e-06 ***
## Fish_size     0.6622     0.2796   2.369   0.0179 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Fish_size -0.926

r2(my_glmC24_FL1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.226
##      Marginal R2: 0.108

5e): Model selection: Squidpops 2h

Temperature

my_glmS2_T1<- glmer(Consumed_individuals~Temperature + (1|Site), data=Squidpops2_stand, family = poisson)

check_overdispersion(my_glmS2_T1)

## # Overdispersion test
## 
##        dispersion ratio =   1.054
##   Pearson's Chi-Squared = 123.348
##                 p-value =   0.326

## No overdispersion detected.

check_model(my_glmS2_T1)

AICc(my_glmS2_T1)

## [1] 460.4034

summary(my_glmS2_T1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.2    468.6   -227.1    454.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7049 -0.8189 -0.1183  0.6789  2.3417 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5106   0.7145  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.4209     0.5139   0.819    0.413
## Temperature   0.3832     0.7173   0.534    0.593
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.903

r2(my_glmS2_T1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.612
##      Marginal R2: 0.016

Salinity

my_glmS2_S1<- glmer(Consumed_individuals~Salinity_ppt + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_S1)

AICc(my_glmS2_S1)

## [1] 460.4596

summary(my_glmS2_S1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.3    468.6   -227.1    454.3      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72415 -0.80815 -0.08599  0.69023  2.37161 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5202   0.7212  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.9474     0.6296   1.505    0.132
## Salinity_ppt  -0.4192     0.8834  -0.475    0.635
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.935

r2(my_glmS2_S1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.615
##      Marginal R2: 0.012

Dissolved oxygen

my_glmS2_Ox1<- glmer(Consumed_individuals~D_oxygen_mgL + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_Ox1)

AICc(my_glmS2_Ox1)

## [1] 457.1796

summary(my_glmS2_Ox1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    457.0    465.3   -225.5    451.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72360 -0.79117 -0.07451  0.73700  2.40980 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3718   0.6098  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    0.1137     0.3581   0.318   0.7508  
## D_oxygen_mgL   1.2713     0.6529   1.947   0.0515 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.844

r2(my_glmS2_Ox1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.606
##      Marginal R2: 0.166

Visibility

my_glmS2_V1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_V1)

AICc(my_glmS2_V1)

## [1] 460.1753

summary(my_glmS2_V1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.0    468.3   -227.0    454.0      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7021 -0.8107 -0.1061  0.7043  2.4274 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5109   0.7148  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    0.4194     0.4157   1.009    0.313
## V_Visibility   0.5814     0.8141   0.714    0.475
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.847

r2(my_glmS2_V1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.617
##      Marginal R2: 0.028

Community composition

my_glmS2_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_PCO1)

AICc(my_glmS2_PCO1)

## [1] 457.9922

summary(my_glmS2_PCO1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    457.8    466.1   -225.9    451.8      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.74367 -0.78077 -0.08552  0.71274  2.46497 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.3944   0.628   
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.3370     0.2866   1.176   0.2396  
## PCo1          0.9343     0.5457   1.712   0.0869 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.726

r2(my_glmS2_PCO1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.604
##      Marginal R2: 0.133

Kelp cover

my_glmS2_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_K1)

AICc(my_glmS2_K1)

## [1] 459.4048

summary(my_glmS2_K1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    459.2    467.6   -226.6    453.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6926 -0.8389 -0.1096  0.6556  2.4268 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4541   0.6739  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   0.4384     0.2948   1.487    0.137
## Kelp          0.8480     0.7294   1.163    0.245
## 
## Correlation of Fixed Effects:
##      (Intr)
## Kelp -0.702

r2(my_glmS2_K1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.606
##      Marginal R2: 0.067

Bare rock cover

my_glmS2_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_BR1)

AICc(my_glmS2_BR1)

## [1] 460.4304

summary(my_glmS2_BR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.2    468.6   -227.1    454.2      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7003 -0.8230 -0.1211  0.6796  2.3980 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5128   0.7161  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.7914     0.3284   2.410    0.016 *
## Bare_rock    -0.3255     0.6415  -0.507    0.612  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.738

r2(my_glmS2_BR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.612
##      Marginal R2: 0.014

Understorey algae cover

my_glmS2_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_UA1)

AICc(my_glmS2_UA1)

## [1] 460.6711

summary(my_glmS2_UA1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.5    468.8   -227.2    454.5      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.71045 -0.81805 -0.09991  0.68648  2.36196 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5282   0.7267  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)  
## (Intercept)        0.68441    0.27845   2.458    0.014 *
## Understorey_algae -0.06958    0.62237  -0.112    0.911  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.592

r2(my_glmS2_UA1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.614
##      Marginal R2: 0.001

Species richness

my_glmS2_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_SR1)

AICc(my_glmS2_SR1)

## [1] 460.6386

summary(my_glmS2_SR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.4    468.8   -227.2    454.4      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.71377 -0.81473 -0.09892  0.68827  2.36015 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5293   0.7275  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   0.7440     0.4298   1.731   0.0835 .
## Mean_SR      -0.1641     0.7736  -0.212   0.8320  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.853

r2(my_glmS2_SR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.615
##      Marginal R2: 0.003

Abundance of invertebrate predators

my_glmS2_I1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_I1)

AICc(my_glmS2_I1)

## [1] 460.0393

summary(my_glmS2_I1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    459.8    468.2   -226.9    453.8      117 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7185 -0.8250 -0.1281  0.7024  2.4144 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4876   0.6983  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                        Estimate Std. Error z value Pr(>|z|)  
## (Intercept)              0.9020     0.3548   2.542    0.011 *
## Invertebrate_predators  -0.5317     0.6507  -0.817    0.414  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Invrtbrt_pr -0.792

r2(my_glmS2_I1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.609
##      Marginal R2: 0.036

Abundance of fish predators

my_glmS2_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_PF1)

AICc(my_glmS2_PF1)

## [1] 460.5443

summary(my_glmS2_PF1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    460.3    468.7   -227.2    454.3      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.70479 -0.82195 -0.09594  0.67244  2.34183 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.5245   0.7242  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      0.5755     0.3320   1.733    0.083 .
## Predatory_fish   0.2390     0.6392   0.374    0.709  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Predtry_fsh -0.739

r2(my_glmS2_PF1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.615
##      Marginal R2: 0.008

Fish length

my_glmS2_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Squidpops2_stand, family = poisson)

check_model(my_glmS2_FL1)

AICc(my_glmS2_FL1)

## [1] 458.2359

summary(my_glmS2_FL1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
##    Data: Squidpops2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    458.0    466.4   -226.0    452.0      117 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72996 -0.80127 -0.09444  0.71113  2.33612 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.4281   0.6543  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   1.4069     0.4953   2.841   0.0045 **
## Fish_size    -1.2026     0.7565  -1.590   0.1119   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Fish_size -0.911

r2(my_glmS2_FL1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.617
##      Marginal R2: 0.124

5f) Model selection: Squidpops 24h

Temperature

my_glmS24_T1<- glmer(Consumed_individuals~Temperature+ (1|Site), data=Squidpops24_stand, family = poisson)

## boundary (singular) fit: see help('isSingular')

my_glmS24_T1_v2<- glm(Consumed_individuals~Temperature, data=Squidpops24_stand, family = poisson)

check_overdispersion(my_glmS24_T1_v2)

## # Overdispersion test
## 
##        dispersion ratio =  0.331
##   Pearson's Chi-Squared = 39.057
##                 p-value =      1

## No overdispersion detected.

check_model(my_glmS24_T1_v2)

AICc(my_glmS24_T1_v2)

## [1] 443.5292

summary(my_glmS24_T1_v2)

## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1227     0.1166   9.632  < 2e-16 ***
## Temperature   0.4879     0.1575   3.098  0.00195 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 51.673  on 118  degrees of freedom
## AIC: 443.43
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_T1_v2)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.200

Salinity

my_glmS24_S1<- glm(Consumed_individuals~Salinity_ppt, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_S1)

AICc(my_glmS24_S1)

## [1] 447.5637

summary(my_glmS24_S1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Salinity_ppt, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.1271     0.1434   7.857 3.93e-15 ***
## Salinity_ppt   0.4670     0.1956   2.387    0.017 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 55.707  on 118  degrees of freedom
## AIC: 447.46
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_S1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.122

Dissolved oxygen

my_glmS24_O1<- glm(Consumed_individuals~D_oxygen_mgL, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_O1)

AICc(my_glmS24_O1)

## [1] 452.8848

summary(my_glmS24_O1)

## 
## Call:
## glm(formula = Consumed_individuals ~ D_oxygen_mgL, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.39170    0.08033  17.324   <2e-16 ***
## D_oxygen_mgL  0.12467    0.14949   0.834    0.404    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.028  on 118  degrees of freedom
## AIC: 452.78
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_O1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.014

Visibility

my_glmS24_V1<- glm(Consumed_individuals~V_Visibility, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_V1)

AICc(my_glmS24_V1)

## [1] 452.1049

summary(my_glmS24_V1)

## 
## Call:
## glm(formula = Consumed_individuals ~ V_Visibility, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.53164    0.08209  18.659   <2e-16 ***
## V_Visibility -0.20381    0.16916  -1.205    0.228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 60.248  on 118  degrees of freedom
## AIC: 452
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_V1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.030

Community composition

my_glmS24_PC1<- glm(Consumed_individuals~PCo1, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_PC1)

AICc(my_glmS24_PC1)

## [1] 451.8008

summary(my_glmS24_PC1)

## 
## Call:
## glm(formula = Consumed_individuals ~ PCo1, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.50543    0.06171  24.394   <2e-16 ***
## PCo1        -0.16872    0.12776  -1.321    0.187    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.944  on 118  degrees of freedom
## AIC: 451.7
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_PC1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.037

Kelp cover

my_glmS24_K1<- glm(Consumed_individuals~Kelp, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_K1)

AICc(my_glmS24_K1)

## [1] 450.7603

summary(my_glmS24_K1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Kelp, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.37365    0.06285  21.856   <2e-16 ***
## Kelp         0.25829    0.15153   1.705   0.0883 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 58.904  on 118  degrees of freedom
## AIC: 450.66
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_K1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.058

Bare rock cover

my_glmS24_BR1<- glm(Consumed_individuals~Bare_rock, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_BR1)

AICc(my_glmS24_BR1)

## [1] 451.3986

summary(my_glmS24_BR1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Bare_rock, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.51985    0.06556  23.183   <2e-16 ***
## Bare_rock   -0.19684    0.13441  -1.464    0.143    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.542  on 118  degrees of freedom
## AIC: 451.3
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_BR1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.045

Understorey_algae cover

my_glmS24_UA1<- glm(Consumed_individuals~Understorey_algae, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_UA1)

AICc(my_glmS24_UA1)

## [1] 453.5612

summary(my_glmS24_UA1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Understorey_algae, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        1.44308    0.05495  26.261   <2e-16 ***
## Understorey_algae  0.01475    0.12466   0.118    0.906    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.705  on 118  degrees of freedom
## AIC: 453.46
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_UA1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.000

Species richness

my_glmS24_SR1<- glm(Consumed_individuals~Mean_SR, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_SR1)

AICc(my_glmS24_SR1)

## [1] 453.2822

summary( my_glmS24_SR1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Mean_SR, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.40627    0.08759  16.054   <2e-16 ***
## Mean_SR      0.08470    0.15635   0.542    0.588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.426  on 118  degrees of freedom
## AIC: 453.18
## 
## Number of Fisher Scoring iterations: 4

r2( my_glmS24_SR1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.006

Abundance of invertebrate predators

my_glmS24_I1<- glm(Consumed_individuals~Invertebrate_predators, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_I1)

AICc(my_glmS24_I1)

## [1] 453.4745

summary(my_glmS24_I1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Invertebrate_predators, 
##     family = poisson, data = Squidpops24_stand)
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             1.46575    0.07385  19.848   <2e-16 ***
## Invertebrate_predators -0.04313    0.13606  -0.317    0.751    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 61.618  on 118  degrees of freedom
## AIC: 453.37
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_I1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.002

Abundance of predatory fish

my_glmS24_PF1<- glm(Consumed_individuals~Predatory_fish, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_PF1)

AICc(my_glmS24_PF1)

## [1] 451.6715

summary(my_glmS24_PF1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Predatory_fish, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.37803    0.06749  20.418   <2e-16 ***
## Predatory_fish  0.17727    0.12762   1.389    0.165    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.815  on 118  degrees of freedom
## AIC: 451.57
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_PF1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.039

Fish length (cm)

my_glmS24_FL1<- glm(Consumed_individuals~Fish_size, data=Squidpops24_stand, family = poisson)

check_model(my_glmS24_FL1)

AICc(my_glmS24_FL1)

## [1] 451.6657

summary(my_glmS24_FL1)

## 
## Call:
## glm(formula = Consumed_individuals ~ Fish_size, family = poisson, 
##     data = Squidpops24_stand)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.3041     0.1144  11.400   <2e-16 ***
## Fish_size     0.2283     0.1665   1.371     0.17    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 59.809  on 118  degrees of freedom
## AIC: 451.56
## 
## Number of Fisher Scoring iterations: 4

r2(my_glmS24_FL1)

## # R2 for Generalized Linear Regression
##   Nagelkerke's R2: 0.039

Test interaction between temperature and dissolved oxygen

my_glm<- glm(Consumed_individuals~ Temperature * D_oxygen_mgL, data=Squidpops24_stand, family = poisson)
summary(my_glm)

## 
## Call:
## glm(formula = Consumed_individuals ~ Temperature * D_oxygen_mgL, 
##     family = poisson, data = Squidpops24_stand)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.2170     0.1628   7.474 7.78e-14 ***
## Temperature                0.4416     0.2252   1.961   0.0499 *  
## D_oxygen_mgL              -0.4196     0.4782  -0.877   0.3803    
## Temperature:D_oxygen_mgL   0.3665     0.5390   0.680   0.4965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 61.719  on 119  degrees of freedom
## Residual deviance: 50.734  on 116  degrees of freedom
## AIC: 446.49
## 
## Number of Fisher Scoring iterations: 4

AICc(my_glm)

## [1] 446.8355

Part 6: Herbivory intensity: Effect of site specific variables

First, create datasets to be used in the analysis

Herbivory_analysis <- Herbivory_6 %>%
  dplyr::select(-Rope, -Bait, -Camera, -Hours, -Initial_length_mm, -Final_Length_mm, -Scraped, -Survival_proportion)

1) Standardise

Herb_all_num <- Herbivory_analysis %>%
  dplyr::select(-Site, -Bait2,-Consumed_proportion, -Day_length) 
Herb_all_stand <- decostand(Herb_all_num, method="range")
#Test autocorrelation 
corrplot(cor(Herb_all_stand), method = "number")

Herb_all_stand$Site <- Herbivory_analysis$Site
Herb_all_stand$Consumption<- Herbivory_analysis$Consumption
Herb_all_stand$Consumed_proportion<- Herbivory_analysis$Consumed_proportion

6a) Model selection

Temperature

my_glmL_T1<- glmer(Consumed_proportion~  Temperature + (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_overdispersion(my_glmL_T1)

## # Overdispersion test
## 
##  dispersion ratio = 0.876
##           p-value = 0.584

## No overdispersion detected.

check_model(my_glmL_T1)

AICc(my_glmL_T1)

## [1] 444.0932

summary(my_glmL_T1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Temperature + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.1    457.2   -219.0    438.1      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87459 -0.28478 -0.08928  0.12354  3.01044 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.183    2.487   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -0.962      1.740  -0.553     0.58
## Temperature   -2.399      2.515  -0.954     0.34
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temperature -0.883

r2(my_glmL_T1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.671
##      Marginal R2: 0.052

Salinity

my_glmL_S1<- glmer(Consumed_proportion~  Salinity_ppt + (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_S1)

AICc(my_glmL_S1)

## [1] 441.0416

summary(my_glmL_S1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Salinity_ppt + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    441.0    454.1   -217.5    435.0      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.93354 -0.27675 -0.09387  0.13458  2.96049 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 4.546    2.132   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)     1.288      1.846   0.698   0.4853  
## Salinity_ppt   -5.561      2.705  -2.055   0.0398 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Salinty_ppt -0.922

r2(my_glmL_S1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.662
##      Marginal R2: 0.194

Dissolved oxygen

my_glmL_O1<- glmer(Consumed_proportion~  D_oxygen_mgL + (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_O1)

AICc(my_glmL_O1)

## [1] 444.2172

summary(my_glmL_O1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ D_oxygen_mgL + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.2    457.3   -219.1    438.2      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86805 -0.26809 -0.09095  0.11924  2.97058 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.02     2.454   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -3.520      1.451  -2.425   0.0153 *
## D_oxygen_mgL    2.400      2.654   0.904   0.3659  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## D_oxygn_mgL -0.831

r2(my_glmL_O1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.664
##      Marginal R2: 0.050

Visibility

my_glmL_V1<- glmer(Consumed_proportion~ V_Visibility + (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_V1)

AICc(my_glmL_V1)

## [1] 444.6448

summary(my_glmL_V1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ V_Visibility + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.6    457.7   -219.3    438.6      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86293 -0.28102 -0.08826  0.12687  3.02234 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.441    2.538   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    -3.223      1.509  -2.136   0.0327 *
## V_Visibility    1.713      2.851   0.601   0.5480  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## V_Visibilty -0.835

r2(my_glmL_V1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.669
##      Marginal R2: 0.021

Community composition

my_glmL_PC1<- glmer(Consumed_proportion~ PCo1+ (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_PC1)

AICc(my_glmL_PC1)

## [1] 444.2705

summary(my_glmL_PC1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ PCo1 + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.2    457.4   -219.1    438.2      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.88774 -0.28227 -0.08271  0.12286  3.01571 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.24     2.498   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -3.159      1.175  -2.687  0.00721 **
## PCo1           1.863      2.169   0.859  0.39060   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## PCo1 -0.715

r2(my_glmL_PC1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.670
##      Marginal R2: 0.044

Kelp cover

my_glmL_K1<- glmmTMB(Consumed_proportion~Kelp+ (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

check_model(my_glmL_K1)

## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(my_glmL_K1)

## [1] 451.4504

Linear model failed to converge, quadratic model was chosen for comparisons.

summary(my_glmL_K1)

##  Family: binomial  ( logit )
## Formula:          Consumed_proportion ~ Kelp + (1 | Site)
## Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.4    464.5   -222.7    445.4      587 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.716    2.391   
## Number of obs: 590, groups:  Site, 12
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  -2.1796     1.0412  -2.093   0.0363 *
## Kelp         -0.8131     2.5535  -0.318   0.7502  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2(my_glmL_K1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.637
##      Marginal R2: 0.006

Bare rock cover

my_glmL_BR1<- glmer(Consumed_proportion~Bare_rock+ (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_BR1)

AICc(my_glmL_BR1)

## [1] 443.0529

summary(my_glmL_BR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Bare_rock + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    443.0    456.2   -218.5    437.0      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.88882 -0.26670 -0.09483  0.12657  2.98644 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.49     2.343   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   -3.630      1.189  -3.054  0.00226 **
## Bare_rock      3.035      2.136   1.421  0.15545   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## Bare_rock -0.758

r2(my_glmL_BR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.665
##      Marginal R2: 0.107

Understorey algae cover

my_glmL_UA1<- glmer(Consumed_proportion~Understorey_algae+ (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_UA1)

AICc(my_glmL_UA1)

## [1] 437.193

summary(my_glmL_UA1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Understorey_algae + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    437.2    450.3   -215.6    431.2      587 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8504 -0.2511 -0.0863  0.1267  5.0475 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 3.436    1.854   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)   
## (Intercept)        -1.0496     0.7058  -1.487  0.13698   
## Understorey_algae  -6.1317     2.3255  -2.637  0.00837 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Undrstry_lg -0.445

r2(my_glmL_UA1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.714
##      Marginal R2: 0.415

Species richness

my_glmL_SR1<- glmer(Consumed_proportion~Mean_SR+ (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_SR1)

AICc(my_glmL_SR1)

## [1] 444.5391

summary(my_glmL_SR1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ Mean_SR + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    444.5    457.6   -219.2    438.5      587 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.85322 -0.27247 -0.09101  0.12793  2.97373 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 6.295    2.509   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)   -3.364      1.556  -2.162   0.0306 *
## Mean_SR        1.834      2.674   0.686   0.4926  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## Mean_SR -0.849

r2(my_glmL_SR1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.666
##      Marginal R2: 0.027

Abundance of herbivores

my_glmL_H1<- glmmTMB(Consumed_proportion~Herbivores+ (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

check_model(my_glmL_H1)

## `check_outliers()` does not yet support models of class `glmmTMB`.

AICc(my_glmL_H1)

## [1] 451.3434

summary(my_glmL_H1)

##  Family: binomial  ( logit )
## Formula:          Consumed_proportion ~ Herbivores + (1 | Site)
## Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    451.3    464.4   -222.7    445.3      587 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 5.732    2.394   
## Number of obs: 590, groups:  Site, 12
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -2.6513     0.9626  -2.754  0.00588 **
## Herbivores    1.1162     2.4550   0.455  0.64934   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2(my_glmL_H1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.640
##      Marginal R2: 0.012

Phlorotannins

my_glmL_PL1<- glmer(Consumed_proportion~phloro_mean+ (1|Site), data=Herb_all_stand, family = binomial())

## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!

check_model(my_glmL_PL1)

AICc(my_glmL_PL1)

## [1] 431.1052

summary(my_glmL_PL1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Consumed_proportion ~ phloro_mean + (1 | Site)
##    Data: Herb_all_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    431.1    444.2   -212.5    425.1      587 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8432 -0.3082 -0.0483  0.1339 27.0033 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 2.507    1.583   
## Number of obs: 590, groups:  Site, 12
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.08397    0.76929   0.109  0.91309   
## phloro_mean -9.85455    3.39879  -2.899  0.00374 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## phloro_mean -0.689

r2(my_glmL_PL1)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.792
##      Marginal R2: 0.633

Test interaction between temperature and dissolved oxygen

my_glm<- glmer(Consumed_individuals~ Temperature*D_oxygen_mgL + (1|Site), data=Crabs2_stand, family = poisson())
summary(my_glm)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature * D_oxygen_mgL + (1 | Site)
##    Data: Crabs2_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    421.4    435.3   -205.7    411.4      115 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96584 -0.90387  0.06297  0.57760  2.96316 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.2205   0.4696  
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)               -0.6119     0.6037  -1.014    0.311
## Temperature                1.0211     0.8251   1.238    0.216
## D_oxygen_mgL              -0.2223     1.7120  -0.130    0.897
## Temperature:D_oxygen_mgL   1.3775     1.9254   0.715    0.474
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr D_xy_L
## Temperature -0.849              
## D_oxygn_mgL -0.735  0.519       
## Tmprtr:D__L  0.730 -0.674 -0.940

AICc(my_glm)

## [1] 421.9314

my_glm<- glmer(Consumed_individuals~ Temperature*Salinity_ppt + (1|Site), data=Crabs24_stand, family = poisson())

## boundary (singular) fit: see help('isSingular')

summary(my_glm)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature * Salinity_ppt + (1 | Site)
##    Data: Crabs24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##    435.8    449.7   -212.9    425.8      115 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1075 -0.2446  0.1351  0.2772  1.3584 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0        0       
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                0.3196     0.6541   0.489    0.625
## Temperature                0.8991     1.2032   0.747    0.455
## Salinity_ppt               0.8443     1.0544   0.801    0.423
## Temperature:Salinity_ppt  -0.3734     1.7801  -0.210    0.834
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr Slnty_
## Temperature -0.964              
## Salinty_ppt -0.978  0.944       
## Tmprtr:Sln_  0.961 -0.987 -0.972
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

AICc(my_glm)

## [1] 436.2851

my_glm<- glmer(Consumed_individuals~ Temperature*Salinity_ppt + (1|Site), data=Squidpops24_stand, family = poisson())

## boundary (singular) fit: see help('isSingular')

summary(my_glm)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Consumed_individuals ~ Temperature * Salinity_ppt + (1 | Site)
##    Data: Squidpops24_stand
## 
##      AIC      BIC   logLik deviance df.resid 
##      448      462     -219      438      115 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8095 -0.1049  0.1330  0.1700  0.9803 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0        0       
## Number of obs: 120, groups:  Site, 12
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)  
## (Intercept)               0.99054    0.58952   1.680   0.0929 .
## Temperature               0.43957    1.13231   0.388   0.6979  
## Salinity_ppt              0.28427    0.95104   0.299   0.7650  
## Temperature:Salinity_ppt -0.06261    1.66684  -0.038   0.9700  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr Slnty_
## Temperature -0.969              
## Salinty_ppt -0.979  0.948       
## Tmprtr:Sln_  0.966 -0.988 -0.974
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

AICc(my_glm)

## [1] 448.5681

Part 7: Visualisation of best fitted predictors

7a): Crab consumption

Make predictions according to individual models with lowest AICc

Crabs 2 h

p_C2T <- ggpredict(my_glmT1, terms ="Temperature [all]")
p_C2Ox <- ggpredict(my_glmOxC2_1, terms ="D_oxygen_mgL [all]")

Crabs 24 h

p_C24T <- ggpredict(my_glmC24_T1, terms ="Temperature [all]")
p_C24S <- ggpredict(my_glmC24_S1, terms ="Salinity_ppt [all]")
p_C24BR <- ggpredict(my_glmC24_BR1, terms ="Bare_rock [all]")
p_C24K <- ggpredict(my_glmC24_K1, terms ="Kelp [all]")
p_C24FL <- ggpredict(my_glmC24_FL1, terms ="Fish_size [all]")

1) Temperature vs Crabs

ggplot() +
  geom_ribbon(data = p_C2T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "darkslateblue", alpha = .25) +
  geom_line(data = p_C2T, aes(x = x, y =predicted), color =  "darkslateblue", linewidth = 0.6) +
  geom_jitter(data = Crabs2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.008, height = 0.1, size = 1, alpha = 0.5)+
  geom_ribbon(data = p_C24T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24T, aes(x = x, y =predicted), color =  "orange1", linewidth = 0.6) +
  geom_jitter(data = Crabs24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.008, height = 0.1, size = 1, alpha = 0.5)+
  labs(x = "Temperature (°C)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1))

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("TempC.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

2) Salinity vs Crabs

ggplot() +
  geom_ribbon(data = p_C24S, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24S, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  geom_jitter(data = Crabs2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,8), breaks = seq(0, 6, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("SalC.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

3) Oxygen vs Crabs

ggplot() +
  geom_ribbon(data = p_C2Ox, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "darkslateblue", alpha = .25) +
  geom_line(data = p_C2Ox, aes(x = x, y = predicted), color =  "darkslateblue", size = 1) +
  geom_jitter(data = Crabs2_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  geom_jitter(data = Crabs24_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
  labs(x = "Dissolved oxygen (mgL)", y = "Crabs consumed")  + Theme1 +
  scale_x_continuous(breaks = c(0.1082296, 0.384854, 0.66147994, 0.935517247), labels = c("4", "6", "8", "10"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,8), breaks = seq(0, 6, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("OxC.jpeg", units="mm", width=90, height=60, dpi=300)

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

4) Bare rock Crabs 24 h

ggplot() +
  geom_ribbon(data = p_C24BR, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24BR, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Bare rock cover (%)", y = "Consumed crabs")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("BareC.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).

4) Fish length - Crabs 24 h

ggplot() +
  geom_ribbon(data = p_C24FL, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_C24FL, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Crabs24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Fish length (cm)", y = "Consumed crabs")  + Theme1 +
  scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("Fish_length.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

7b): Squidpop consumption

Squidpops 2 h

None

Squidpops 24 h

p_S24T <- ggpredict(my_glmS24_T1, terms ="Temperature [all]")
p_S24S <- ggpredict(my_glmS24_S1, terms ="Salinity_ppt [all]")

1) Temperature vs squidpops

ggplot() +
  geom_ribbon(data = p_S24T , aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_S24T , aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Squidpops24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
   geom_jitter(data = Squidpops2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Temperature (°C)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("TemperatureS.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

2) Salinity vs squidpops

ggplot() +
  geom_ribbon(data = p_S24T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_S24T, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Squidpops24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
   geom_jitter(data = Squidpops2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("SalS.jpeg", units="mm", width=110, height=80, dpi=300)

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_point()`).

Bare rock vs squidpops

ggplot() +
  geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
    geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Bare rock cover (%)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("BRS.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

Fish length vs squidpops

ggplot() +
  geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
    geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
  labs(x = "Fish Length (cm)", y = "Consumed squidpops")  + Theme1 +
  scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
  scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("FLS2.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

7c): Seaweed consumption

p_Lphloro <- ggpredict(my_glmL_PL1, terms ="phloro_mean [all]")
p_L24S <- ggpredict(my_glmL_S1, terms ="Salinity_ppt [all]")
p_L24UA <- ggpredict(my_glmL_UA1, terms ="Understorey_algae [all]")

Phlorotannins vs herbivory

ggplot() +
  geom_ribbon(data = p_Lphloro, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_Lphloro, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= phloro_mean, y=Consumed_proportion), color = "orange1", width=0.02, height = 0.01, size = 1.5, alpha = 0.5)+
  labs(x = "Phlorotannin content (mg"~g^-1*" DW)", y = "Blade consumed proportion")  + Theme1 +
  scale_x_continuous(breaks = c(0.04788192, 0.2840443, 0.5202066, 0.756369, 0.9925314), labels = c("8", "16", "24", "32", "40"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 230 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("PhloroL.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 229 rows containing missing values or values outside the scale range
## (`geom_point()`).

Salinity vs herbivory

ggplot() +
  geom_ribbon(data = p_L24S , aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_L24S, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= Salinity_ppt, y=Consumed_proportion), color = "orange1", width=0.02, height = 0.01, size = 1.5, alpha = 0.5)+
  labs(x = "Salinity (ppt)", y = "Blade consumed proportion")  + Theme1 +
  scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 268 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("SalL.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 240 rows containing missing values or values outside the scale range
## (`geom_point()`).

UA vs herbivory

ggplot() +
  geom_ribbon(data = p_L24UA, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
  geom_line(data = p_L24UA, aes(x = x, y = predicted), color =  "orange1", size = 1) +
  geom_jitter(data = Herb_all_stand, aes(x= Understorey_algae, y=Consumed_proportion), color = "orange1", width=0.03, height = 0.03, size = 1.5, alpha = 0.5)+
  labs(x = "UA cover (%)", y = "Predicted consumed proportion ")  + Theme1 +
  scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
  scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) + 
  theme(legend.position = c(0.1, 0.85))

## Warning: Removed 269 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("UAL.jpeg", units="mm", width=90, height=70, dpi=300)

## Warning: Removed 275 rows containing missing values or values outside the scale range
## (`geom_point()`).

Part 8: Validate sea water parameters

Comparison between data obtained in-situ and long-term data. Sea surface temperature (SST) and dissolved oxygen daily average were extracted from E.U. Copernicus Marine Service Information (doi =10.48670/moi-00169, doi =10.48670/moi-00015, respectively). Salinity data were extracted using the NASA Earth Data Search (doi = SMP20-4U7CS). SST and salinity data were collected between 1 December 2018 and 31 March 2019, which matches the period of the study. Dissolved oxygen data were available for the time period 1 November 2021 and 28 February 2025.

8a) Temperature

Temp <- ggplot() +
  geom_jitter(data = temp_data, aes(x= Latitude*-1, y=analysed_sst), width = 0.08, height = 0.03, size = 0.8, alpha = 0.3) +
  geom_boxplot(data = temp_data, aes(x= Latitude*-1, y=analysed_sst, group =  cut_width(Latitude,1)), alpha = 0.2) +
  geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$Temperature), size=2, color = "red", pch =15)+
  Theme1 +
labs( y="SST (°C)", x="Latitude (°S)") +
  theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x =     element_text(size = 11))

Temp

ggsave("TVar.jpeg", units="mm", width=130, height=100, dpi=300)

8b) Salinity

Sal <- ggplot() +
  geom_jitter(data = Sal_data, aes(x= lat*-1, y=Salinity), width = 0.08, height = 0.03, size = 0.8, alpha = 0.6) +
  geom_boxplot(data = Sal_data, aes(x= lat*-1, y=Salinity, group =  cut_width(lat,1)), alpha = 0.2) +
  geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$Salinity_ppt), size=2, color = "red", pch =15)+
  Theme1 +
labs( y="Salinity (ppt)", x="Latitude (°S)") +
  theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x =     element_text(size = 11))

Sal

ggsave("SalVar.jpeg", units="mm", width=130, height=100, dpi=300)

8c) Oxygen

Transform to mg/L

library(respR)

## Warning: package 'respR' was built under R version 4.4.2

## Registered S3 method overwritten by 'respR':
##   method        from
##   print.inspect pryr

## 
## Attaching package: 'respR'

## The following object is masked from 'package:MASS':
## 
##     select

o2_data$conv <- convert_DO(o2_data$o2, # data to convert
                   from = "umol/L",     # oxygen unit to convert from
                   to = "mg/L",    # oxygen unit to convert to
                   t = 15,            # in C
                   S = 34,            # in ppt
                   P = 1.013)         # in bar

Make plot

ox <- ggplot() +
  geom_jitter(data = o2_data, aes(x= Latitude*-1, y=conv), width = 0.08, height = 0.03, size = 0.8, alpha = 0.2) +
  geom_boxplot(data = o2_data, aes(x= Latitude*-1, y=conv, group =  cut_width(Latitude,1)), alpha = 0.2) +
  geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$D_oxygen_mgL), size=2, color = "red", pch =15)+
  Theme1 +
labs( y="Dissolved oxygen (mg"~ L^-1*")", x="Latitude (°S)") +
  theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x =     element_text(size = 11))

ox

ggsave("o2Var.jpeg", units="mm", width=130, height=100, dpi=300)

Part 9: Principal component analysis (PCA) of predictor variables

In this section, a PCA with all the predictor variables was carried out and contracted with latitude.

Get data ready

Variables <- Phlorotannins_sum  %>%
    left_join(Environmental_sum, by = "Site") %>%
  left_join(Group_abundances_mean, by = "Site") %>%
  left_join(Richness_meanpred, by = "Site") %>%
   left_join( pcoa_df1, by = "Site") %>%
   left_join(Size_mean_wide, by = "Site") %>%
   left_join(Benthic_sum_wide_var, by = "Site")

Variables1 <- dplyr::select(Variables, -Site, -Latitude, -Day_length)

Perform PCA

# Perform PCA
pca_result <- prcomp(Variables1, center = TRUE, scale. = TRUE)
# View PCA summary
summary(pca_result)

## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6     PC7
## Standard deviation     1.9056 1.8232 1.5595 1.2683 1.1188 0.80575 0.67715
## Proportion of Variance 0.2594 0.2374 0.1737 0.1149 0.0894 0.04637 0.03275
## Cumulative Proportion  0.2594 0.4968 0.6705 0.7854 0.8748 0.92121 0.95396
##                            PC8     PC9    PC10    PC11     PC12
## Standard deviation     0.59107 0.42934 0.32003 0.09201 4.97e-15
## Proportion of Variance 0.02495 0.01317 0.00732 0.00060 0.00e+00
## Cumulative Proportion  0.97891 0.99208 0.99940 1.00000 1.00e+00

# Access PCA components
pca_result$rotation  # Loadings (eigenvectors)

##                                PC1         PC2         PC3         PC4
## phloro_mean             0.15312183 -0.38036071  0.16987612 -0.40416794
## D_oxygen_mgL            0.04416129 -0.06947408  0.22900240 -0.29578882
## Salinity_ppt            0.32770211 -0.28725268 -0.30041598  0.01204521
## Temperature             0.33959317 -0.19864160  0.25173653  0.07112323
## V_Visibility           -0.26155336 -0.04655336 -0.46754924 -0.10336321
## Herbivores              0.17945120  0.26495508 -0.32071365 -0.30175664
## Invertebrate_predators  0.23631635  0.35158153  0.08199631 -0.22146569
## Predatory_fish          0.33348022  0.19557271  0.12844257  0.42742255
## Mean_SR                 0.24588446  0.13891112 -0.46658761 -0.10235712
## PCo1                   -0.45336621 -0.12321574  0.20140955  0.11357894
## Fish_size               0.40286210 -0.15837393  0.19130840 -0.19017485
## Kelp                    0.20344908 -0.25734537 -0.07308674  0.49797138
## Bare_rock              -0.02538404  0.41986093  0.34602767 -0.18266636
## Understorey_algae      -0.11785619 -0.43843959 -0.02128394 -0.26379807
##                                PC5          PC6         PC7         PC8
## phloro_mean            -0.22928655 -0.061934912  0.16030316  0.06977640
## D_oxygen_mgL            0.71917764 -0.141970506  0.24870627  0.01063205
## Salinity_ppt           -0.04621847  0.088594787 -0.25384253  0.37791097
## Temperature             0.34233239  0.366339285  0.01239502  0.29709948
## V_Visibility            0.28277131 -0.265296643 -0.32346405 -0.03883120
## Herbivores             -0.18909046 -0.075864110  0.65389541  0.11582022
## Invertebrate_predators  0.07850895 -0.516299347 -0.21984644  0.32372054
## Predatory_fish         -0.04690808 -0.353793051 -0.09874146  0.21766508
## Mean_SR                 0.31766340  0.122121426 -0.05131298 -0.29528597
## PCo1                    0.18079042 -0.201314998  0.07870739  0.07955889
## Fish_size              -0.07309108 -0.101363639 -0.28571931 -0.59602278
## Kelp                    0.08064145 -0.390533765  0.39301185 -0.31521002
## Bare_rock              -0.07793862 -0.008588723 -0.09287449 -0.20915738
## Understorey_algae      -0.18785568 -0.384416932 -0.06340740  0.07047988
##                                PC9         PC10         PC11         PC12
## phloro_mean            -0.19914246 -0.302226749  0.547906843 -0.325191416
## D_oxygen_mgL            0.28917754 -0.203944032 -0.080842398 -0.005764669
## Salinity_ppt            0.29382455 -0.354717840  0.055501156  0.535005519
## Temperature            -0.10134383  0.454952128 -0.004250124 -0.076982913
## V_Visibility            0.07492272 -0.003118584  0.171242875 -0.276713426
## Herbivores              0.13103355 -0.015917972 -0.278008088  0.010513205
## Invertebrate_predators -0.52190028  0.112816781  0.002910929  0.218806222
## Predatory_fish          0.38672379 -0.153625509 -0.023726204 -0.490900382
## Mean_SR                -0.06100320  0.112910544  0.207166319 -0.035055321
## PCo1                   -0.11580691 -0.293793915 -0.143554874  0.226404458
## Fish_size              -0.04289802 -0.178281745 -0.435374536  0.012026704
## Kelp                   -0.11866981  0.127147103  0.290625846  0.278306208
## Bare_rock               0.44193909  0.139799385  0.473690956  0.325001830
## Understorey_algae       0.32715632  0.575421427 -0.145109455  0.028404422

pca_result$x         # Principal components (scores)

##               PC1          PC2         PC3        PC4         PC5        PC6
##  [1,]  1.01934232 -2.465942186 -0.05856192 -1.1805419 -0.74403145 -0.6008792
##  [2,]  1.64364575 -0.247652503  0.33059729  2.8574056  0.41942831 -0.9032240
##  [3,]  0.65673090 -1.475702438 -0.83220277  0.1380973 -0.43291297  1.0595816
##  [4,]  0.10429367  3.053015250 -2.55270043 -1.8166224  0.04174398 -0.5924529
##  [5,]  0.99165334  0.069541318  0.37106187 -0.4478629 -1.35364440  0.5912784
##  [6,] -4.14190953 -0.901118939 -2.06185691  1.0229257 -0.76843221 -0.5500265
##  [7,]  1.57351825  1.793784930  1.24878295 -0.1266223 -0.94011881  0.5183558
##  [8,]  0.52079420  2.322458024  1.45717978  0.2800948  0.28004323 -1.0837141
##  [9,] -0.02668262 -0.009488045 -0.64428859 -0.2340490  2.98842103  0.7665162
## [10,] -0.39575077 -2.747177050  1.78174925 -1.6003581  0.71592698 -0.7627492
## [11,] -3.50506614  1.243557346  2.39160564  0.2027377 -0.10405266  0.8614225
## [12,]  1.55943063 -0.635275706 -1.43136615  0.9047954 -0.10237103  0.6958914
##              PC7         PC8         PC9        PC10        PC11          PC12
##  [1,]  0.2498895 -0.72786447  0.49737574  0.47748733  0.08635639 -4.642120e-15
##  [2,]  0.5267490 -0.03407682 -0.56423696  0.21510138  0.01157954 -5.436623e-15
##  [3,] -0.8074396  0.10588014 -0.19258732  0.36593723 -0.18791464 -5.256212e-15
##  [4,]  0.6350624 -0.05283392 -0.36977732  0.09519209 -0.04917014 -4.992534e-15
##  [5,] -0.6701076 -0.56860844 -0.67721358 -0.40815778  0.11667248 -4.359360e-15
##  [6,] -0.4861734  0.48875365  0.14721247 -0.08360601  0.05770054 -5.280498e-15
##  [7,]  0.1000942  1.32580148  0.29863352  0.19812510  0.07925996 -4.798245e-15
##  [8,] -1.1046837 -0.41762546  0.51534582 -0.17420288 -0.07669713 -4.576201e-15
##  [9,] -0.4462995 -0.02123266  0.04596192  0.08061882  0.10601955 -4.909267e-15
## [10,]  0.3077744  0.66291487 -0.28732577 -0.34383257 -0.06506740 -4.635181e-15
## [11,]  0.8137796 -0.57304152 -0.01998176  0.14372574 -0.03156541 -4.538037e-15
## [12,]  0.8813548 -0.18806684  0.60659323 -0.56638845 -0.04717373 -4.943962e-15

# Scree plot
fviz_eig(pca_result, addlabels = TRUE)

Obtain variables for ggplot

pca_scores <- as.data.frame(pca_result$x)
pca_scores$Site <- Variables$Site
pca_scores$Latitude <- Variables$Latitude
pca_loadings <- as.data.frame(pca_result$rotation)
pca_loadings$Variable <- rownames(pca_loadings)



# Rename variable labels in the pca_loadings data frame
pca_loadings$Variable <- dplyr::recode(pca_loadings$Variable,
  V_Visibility = "Visibility",
  Understorey_algae = "UA cover",
  Bare_rock = "Bare rock cover",
  D_oxygen_mgL = "Dissolved oxygen",
  phloro_mean = "Phlorotannins",
  Kelp = "Kelp cover",
  PCo1 = "PCo1",
  Temperature = "Temperature",
  Fish_size = "Fish length",
  Salinity_ppt = "Salinity",
  Predatory_fish = "Predatory fish",
  Invertebrate_predators = "Invertebrate predators",
  Mean_SR = "Species richness",
  Herbivores = "Herbivores"
)

Visualise

ggplot() +
  # Add PCA scores (samples) with latitude gradient
  geom_point(data = pca_scores, aes(x = PC1, y = PC2, color = Latitude*-1), size = 3) +
  scale_color_viridis_c(option = "plasma", trans = "reverse") +
  # Add PCA loadings (variables as arrows)
  geom_segment(data = pca_loadings,
               aes(x = 0, y = 0, xend = PC1 * 5, yend = PC2 * 5), 
               arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  # Add variable labels
  geom_text(data = pca_loadings, 
            aes(x = PC1 * 5, y = PC2 * 5, label = Variable), 
            color = "black", vjust = -0.5) +
  # Customize plot appearance
  labs(      x = "Principal Component 1",
       y = "Principal Component 2",
       color = "Latitude (°S)") +
  Theme1

ggsave("PCA.jpeg", units="mm", width=130, height=120, dpi=300)

Evaluate statistical effects of latitude

# Load vegan package
library(vegan)

# Create a distance matrix from PCA scores
pca_dist <- dist(pca_scores[, c("PC1", "PC2")])


# Perform PERMANOVA
permanova_result <- adonis2(pca_dist ~ Latitude, data = pca_scores, permutations = 999)
print(permanova_result)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = pca_dist ~ Latitude, data = pca_scores, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)    
## Latitude  1   33.467 0.43742 7.7754  0.001 ***
## Residual 10   43.042 0.56258                  
## Total    11   76.510 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Part 10: Size estimations summary table

A summary table of the size estimations by species and site was created.

Size_summary <- Size_metadata %>%
  group_by(Species) %>%
  summarise (N_individuals = n(), mean_Size_sp = mean(Size_cm), sd = sd(Size_cm), Max_size = max(Size_cm), Min_size = min(Size_cm))

kable(Size_summary)

Species	N_individuals	mean_Size_sp	sd	Max_size	Min_size
Acanthistius_pictus	3	16.000000	1.0000000	17	15
Anisotremus_scapularis	4	21.250000	2.9860788	25	18
Aplodactylus_punctatus	54	27.944444	6.5196848	40	14
Auchenionchus_critinus	1	12.000000	NA	12	12
Auchenionchus_sp.	14	19.071429	6.0949448	32	11
Auchenionchus_variolosus	3	28.666667	5.6862407	35	24
Bovichtus_chilensis	3	17.666667	9.0184995	27	9
Calliclinus_geniguttatus	19	10.157895	2.3632172	14	4
Cheilodactylus_variegatus	125	22.560000	7.0715240	38	6
Cheilotrema_fasciatum	1	28.000000	NA	28	28
Chromis_crusma	85	13.611765	4.2457096	26	5
Congiopodus_peruvianus	2	13.500000	0.7071068	14	13
Girella_laevifrons	11	26.272727	7.3087743	38	15
Graus_nigra	15	25.066667	4.9923751	37	19
Hemilutjanus_macrophthalmos	2	20.000000	2.8284271	22	18
Hepatus_chilensis	1	6.000000	NA	6	6
Homalaspis_plana	20	9.800000	1.9084301	12	4
Hypsoblennius_sordidus	5	7.800000	1.0954451	9	6
Isacia_conceptionis	7	25.571429	2.5071327	30	22
Labrisomus_philippii	18	28.055556	6.9745101	38	15
Nexilosus_latifrons	4	23.000000	5.7154761	30	17
Paralabrax_humeralis	14	26.571429	8.4826908	44	15
Paralichthys_microps	4	32.500000	23.1012265	65	16
Paraxanthus_barbiger	12	6.250000	1.6583124	9	4
Patagonotothen_brevicauda	10	8.100000	1.3703203	10	6
Pinguipes_chilensis	81	28.370370	8.5197483	51	13
Prolatilus_jugularis	20	24.000000	6.9357958	39	10
Romaleon_setosum	17	8.058824	2.2768012	12	5
Scartichthys_spp.	176	13.670454	5.0927884	34	4
Schroederichthys_chilensis	3	51.666667	2.8867513	55	50
Sebastes_oculatus	12	12.333333	7.1137937	24	4
Semicossyphus_darwini	6	26.500000	6.4420494	36	20
Taliepus_sp.	27	9.444444	1.8045526	12	5

write.csv(Size_summary, "Size_summary2.csv")

Part 11: R Session Information

In this part, information about the operating system and R packages attached during the production of this document can be found

sessionInfo() %>% pander

R version 4.4.1 (2024-06-14 ucrt)

Platform: x86_64-w64-mingw32/x64

locale: LC_COLLATE=English_Australia.utf8, LC_CTYPE=English_Australia.utf8, LC_MONETARY=English_Australia.utf8, LC_NUMERIC=C and LC_TIME=English_Australia.utf8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: respR(v.2.3.3), fitdistrplus(v.1.2-2), survival(v.3.6-4), factoextra(v.1.0.7), lme4(v.1.1-35.4), Matrix(v.1.7-0), glmmTMB(v.1.1.10), mgcv(v.1.9-1), nlme(v.3.1-167), geosphere(v.1.5-18), MuMIn(v.1.48.4), MASS(v.7.3-60.2), corrplot(v.0.95), car(v.3.1-2), carData(v.3.0-5), knitr(v.1.49), performance(v.0.13.0), ggeffects(v.1.7.0), pander(v.0.6.5), patchwork(v.1.3.0), vegan(v.2.6-6.1), lattice(v.0.22-6), permute(v.0.9-7), summarytools(v.1.0.1), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): rstudioapi(v.0.16.0), DHARMa(v.0.4.6), jsonlite(v.1.9.0), shape(v.1.4.6.1), datawizard(v.1.0.0), magrittr(v.2.0.3), TH.data(v.1.1-2), estimability(v.1.5.1), magick(v.2.8.3), farver(v.2.1.2), nloptr(v.2.1.0), rmarkdown(v.2.29), ragg(v.1.3.2), vctrs(v.0.6.5), minqa(v.1.2.7), base64enc(v.0.1-3), rstatix(v.0.7.2), htmltools(v.0.5.8.1), haven(v.2.5.4), broom(v.1.0.6), marelac(v.2.1.11), sass(v.0.4.9), bslib(v.0.9.0), plyr(v.1.8.9), sandwich(v.3.1-1), emmeans(v.1.10.2), zoo(v.1.8-12), cachem(v.1.1.0), TMB(v.1.9.16), lifecycle(v.1.0.4), pkgconfig(v.2.0.3), sjlabelled(v.1.2.0), R6(v.2.6.1), fastmap(v.1.2.0), rbibutils(v.2.2.16), digest(v.0.6.35), numDeriv(v.2016.8-1.1), colorspace(v.2.1-0), textshaping(v.0.4.0), ggpubr(v.0.6.0), labeling(v.0.4.3), timechange(v.0.3.0), abind(v.1.4-5), compiler(v.4.4.1), withr(v.3.0.2), backports(v.1.5.0), ggsignif(v.0.6.4), gsw(v.1.2-0), roll(v.1.1.7), tools(v.4.4.1), glue(v.1.7.0), grid(v.4.4.1), checkmate(v.2.3.1), cluster(v.2.1.6), reshape2(v.1.4.4), see(v.0.10.0), generics(v.0.1.3), gtable(v.0.3.6), tzdb(v.0.4.0), data.table(v.1.15.4), hms(v.1.1.3), sp(v.2.1-4), xml2(v.1.3.6), ggrepel(v.0.9.6), pillar(v.1.10.1), splines(v.4.4.1), pryr(v.0.1.6), oce(v.1.8-3), tidyselect(v.1.2.1), reformulas(v.0.4.0), stats4(v.4.4.1), xfun(v.0.51), rapportools(v.1.1), matrixStats(v.1.5.0), stringi(v.1.8.4), yaml(v.2.3.8), boot(v.1.3-30), evaluate(v.1.0.3), codetools(v.0.2-20), tcltk(v.4.4.1), seacarb(v.3.3.3), SolveSAPHE(v.2.1.0), cli(v.3.6.2), RcppParallel(v.5.1.10), xtable(v.1.8-4), systemfonts(v.1.1.0), Rdpack(v.2.6), munsell(v.0.5.1), jquerylib(v.0.1.4), Rcpp(v.1.0.12), coda(v.0.19-4.1), parallel(v.4.4.1), bayestestR(v.0.15.1), viridisLite(v.0.4.2), mvtnorm(v.1.2-5), scales(v.1.3.0), insight(v.1.0.2), crayon(v.1.5.3), rlang(v.1.1.4) and multcomp(v.1.4-26)

Environmental and local habitat variables as predictors of trophic interactions in subtidal rocky reefs along the SE Pacific coast

12-08-2024

Clear the environment

Load packages

Part 1: Importing and examining data

1a) Import data

1b) Inspect structure of data frames

Part 2: Manipulating datasets

2a) Environmental data

2b): Benthic data

2c): Fish length

2d) Species list (abundances and richness)

2e): Community composition

2f): Phlorotannins

Part 3: Latitudinal patterns in explanatory variables

3a): Explore latitudinal relationships

1) Temperature vs latitude

2) Salinity vs latitude

3) Oxygen vs latitude

4) Visibility vs latitude

5) PCOA vs latitude

6) Kelp vs latitude

7) Bare_rock vs latitude

8) Understorey_algae vs latitude

9) Predatory fish abundance vs latitude

10) Herbivorous fish abundance vs latitude

11) Invertebrate predators abundance vs latitude

12) Species richness vs latitude

13) Fish length vs latitude

14) Phlorotannins vs latitude

Final plots

Part 4: Latitudinal patterns in consumption

4a): Predation vs latitude

1) Crabs 2h

2) Crabs 24h

3) Crabs 2 and 24 h prediction plot

4) Squidpops 2h

5) Squidpops 24h

3) Squidpops 2 and 24 h prediction plot

4b): Herbivory vs latitude

Part 5: Predation intensity: Effect of site-specific variables

5a): Standardisation and test autocorrelation

5b) Model selection: Crabs 2h

5c) Model selection: Crabs 24h

5e): Model selection: Squidpops 2h

5f) Model selection: Squidpops 24h

Part 6: Herbivory intensity: Effect of site specific variables

6a) Model selection

Part 7: Visualisation of best fitted predictors

7a): Crab consumption

7b): Squidpop consumption

7c): Seaweed consumption

Part 8: Validate sea water parameters

8a) Temperature

8b) Salinity

8c) Oxygen

Part 9: Principal component analysis (PCA) of predictor variables

Part 10: Size estimations summary table

Part 11: R Session Information